A whole-genome transcriptome analysis of articular chondrocytes in secondary osteoarthritis of the hip

Objective To date, exhaustive gene expression analyses of chondrocytes in hip osteoarthritis (OA) have yielded specific gene expression patterns. No study has reported on the exhaustive transcriptome of secondary hip OA based on acetabular dysplasia in a Japanese population, while previous reports have focused on primary or idiopathic hip OA in Caucasian populations. This study aims to search for specific gene expression patterns of secondary hip OA chondrocytes by transcriptome analysis. Design Human articular cartilage was obtained from femoral heads following hemiarthroplasty for femoral neck fracture (N = 8; non-OA) and total hip arthroplasty for secondary hip OA (N = 12). Total RNA was extracted from the articular cartilage and submitted for microarray analysis. The obtained data were used to perform gene expression analysis, GO enrichment analysis and pathway analysis and were compared with data from primary hip OA in Caucasian populations in the literature. Results We identified 888 upregulated (fold change: FC ≥ 2) and 732 downregulated (FC ≤ 0.5) genes in hip OA versus non-OA chondrocytes, respectively. Only 10% of upregulated genes were common between the secondary and primary OA. The newly found genes prominently overexpressed in the secondary hip OA chondrocytes were DPT, IGFBP7, and KLF2. Pathway analysis revealed extracellular matrix (ECM)-receptor interaction as an OA-related pathway, which was similar to previous reports in primary hip OA. Conclusions This is the first study to report the genome-wide transcriptome of secondary hip OA chondrocytes and demonstrates new potential OA-related genes. Gene expression patterns were different between secondary and primary hip OA, although the results of pathway and functional analysis were similar.


Introduction
Osteoarthritis (OA) is a joint disease characterized by irreversible, degenerative changes in joint components such as articular cartilage, synovium, and subchondral bone [1]. It is one of the leading causes of disability occurring mainly in the elderly due to severe joint pain and deformity especially in weight-bearing joints such as the knee and hip. Approximately 10-15% of the world's population suffers from symptomatic OA [2]. In Japan, 25 and 12 million people were estimated to suffer from knee OA and hip OA, respectively [3,4]. The prevalence of OA increases yearly as the number of elderly rise. Although OA is a global "common disease" as mentioned above, no definitive therapy has been developed to alleviate degraded joint components. Once joint destruction is at the end stage, surgeries such as arthroplasty or joint replacement are the only solution to relieve joint pain and deformity. These days, medical services for OA are giving a huge socioeconomical impact on increasing medical costs [5].
OA develops due to various causes including aging, obesity, mechanical load and genetic factors. Previous studies have investigated a genetic influence in OA [6] and suggested that genetic susceptibility of OA varies by joints, such as the hip or knee [7]. Furthermore, hip OA is more influenced by genetic factors than are other joints. Recently, genetic research has developed rapidly with novel sequencing technologies and microarray analysis [8]. Exhaustive transcriptome analysis of chondrocytes in primary hip OA has uncovered some specific gene expression patterns. Idiopathic OA, also known as primary OA is the dominant form of hip OA in Caucasians [9]. Meanwhile in Japan, approximately 80-90% of hip OA is considered to be secondary following acetabular dysplasia (AD) [10,11]. Considering the mechanisms of onset and progression, primary and secondary hip OA should have different pathomechanisms. In this respect, gene expression patterns in the chondrocytes of primary hip OA in previous reports are assumed to be different from those in secondary OA Japanese patients. Most exhaustive gene analyses have been targeted to primary OA [12,13]; however, no exhaustive gene analysis has been conducted yet in patients with secondary hip OA. Therefore, this study aims to identify specific gene expression patterns with microarray analysis of the articular chondrocytes on secondary hip OA in Japanese population. Furthermore, differences or similarities in the transcriptome of OA chondrocytes between the primary and secondary are discussed and compared with previous reports.

Ethics statement
This study was approved by the Ethical Review Board of Tohoku University Graduate School of Medicine. Informed consent was obtained from all patients for the use of study data.

Human articular cartilage harvesting
Human articular cartilage pieces were obtained after hemiarthroplasty following fracture of the neck of femur (#NOF: non-OA), or total hip arthroplasty for secondary hip OA. Surgeries were performed at affiliated hospitals of the Department of Orthopaedic Surgery, Tohoku University, Sendai, Japan. Clinical findings for each patient were recorded by survey form including medical history, medication history, and radiographic classification of OA [14]. Radiological evaluation of hip OA was performed using the Japanese Orthopaedic Association staging criteria. The center-edge angle, the angle formed by a vertical line and a line drawn from the center of the femoral head to the edge of the acetabulum on plain X-ray, <20 degrees was defined as AD of the hip [10]. All OA patients in this study were in the terminal stage of secondary hip OA based on AD. Patients with primary hip OA, rheumatoid arthritis, or osteonecrosis of the hip were excluded from the study.

RNA extraction
Articular cartilage pieces were harvested from femoral heads within 6 hours of surgery. Non-OA cartilage samples were obtained from the middle to deep zone of microscopically intact, non weight-bearing part of #NOF femoral heads as shown in the previous study [15]. OA cartilage samples were harvested from severely eroded areas of femoral head cartilage next to the weight-bearing area (Fig 1).
Cartilage pieces were finely diced and suspended in QIAzol 1 Lysis Reagent (Qiagen, Crawley, UK). The suspension was homogenized on ice with TissueRuptor™ (Qiagen, Crawley, UK) for efficient disruption of the cartilage pieces. Total RNA was directly extracted from the supernatant prepared above using the Qiagen RNeasy Lipid Tissue Mini Kit 1 (Qiagen, Crawley, UK) according to the manufacturer's instructions. Concentration and quality of total RNA were assessed by NanoDrop™ Lite (Thermo Fisher Scientific, Inc., Waltham, MA, USA) and Agilent 2100 Bioanalyzer 1 (Agilent Technologies, Santa Clara, CA, USA), respectively. Extracted total RNA was also reverse transcribed with High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, CA, USA) immediately following extraction for quantitative real-time polymerase chain reaction (qRT-PCR).

Transcriptome analysis with microarray
Extracted total RNA from 8 non-OA and 12 OA chondrocytes with RNA Integrity Number ! 6.5 underwent microarray analysis using 3D-Gene Human Oligo Chip 25k (Toray Industries Inc., Tokyo, Japan). Demographic data of the samples and quality of RNA are shown in Table 1. Extracted total RNA was labeled with Cy5 using the Amino Allyl MessageAmp™ II aRNA Amplification Kit (Applied Biosystems, CA, USA). The Cy5-labeled aRNA pools were applied to a hybridization buffer and hybridized for 16 h following the supplier's protocols (www.3d-gene.com). Fluorescent signals were scanned using a 3D-Gene Scanner (Toray Industries Inc., Tokyo, Japan) and analyzed with 3D-Gene Extraction software (Toray Industries Inc., Tokyo, Japan). A global normalization method was used to adjust the median of all detected signal intensities to 25. The genes were selected for further analysis when (1) expression with a mean fold change (FC) of 0.5 or ! 2 between the non-OA and OA groups was observed and (2) the round-robin comparison of the gene expression between the 8 non-OA and 12 OA samples (96 combinations in total) showed a significant difference in > 80% (! 77 of 96) of the combinations. The analysis aims to narrow down the genes with higher specificity without affecting the average value of gene expression according to the previous report [16]. A p-value < 0.01 was considered statistically significant.
Gene Ontology (GO) functional enrichment analysis for the differentially expressed genes was performed by Gene Set Enrichment Analysis software (http://software.broadinstitute.org/ gsea/index.jsp) [17]. The gene sets were separated according to the GO terms for biological processes, cellular components, and molecular functions. A p-value < 0.01 and false discovery rate (FDR) q-values 0.25 were used to filter results. Pathway analysis was performed using GeneCodis tools (http://genecodis.cnb.csic.es/) based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database with a p-value < 0.05 considered statistically significant [18].

Quantitative real-time polymerase chain reaction (qRT-PCR)
To replicate the gene expression profiles in the microarray analysis, total RNA extracted from 6 independent non-OA and 6 OA samples were analyzed by qRT-PCR. Twelve differently expressed genes in OA chondrocytes, ASPN, COL1A2, COL2A1, COL3A1, COL5A2, KLF2, MXRA5, OGN, SPARC, TGFBI, and TNFAIP6, were selected according to our microarray data with the following conditions: (1) mean FC ! 10 in OA compared to non-OA groups, (2) expression signals ! 1000 in OA, and (3) > 80% (! 77 of 96 combinations) of the roundrobin comparisons between 8 non-OA and 12 OA samples showed significantly increased expression in OA. Relative quantification of gene expression was performed with an Applied Biosystems 1 StepOnePlus™ Real-time PCR System (Applied Biosystems, CA, USA) using Taq-Man 1 Gene Expression Assays (Applied Biosystems). Reactions were performed in duplicate with glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as an internal control gene. A 20-μL reaction mixture was prepared for each reaction containing 1 μL of complementary DNA, 10 μL of TaqMan 1 Universal PCR Master Mix (Applied Biosystems), 1 μL of premixed gene assay mix with 5 μM of TaqMan 1 probe and 18 μM for each primer, and 8 μL of nuclease-free water. Thermocycler conditions consisted of an initial activation at 95˚C for 10 minutes, followed by a 2-step PCR program at 95˚C for 15 seconds and 60˚C for 60 seconds for 40 cycles. The 2 -ΔΔCt method was used for relative quantification of gene expression. A dissociation curve was obtained after each qPCR run to validate reactions. Data are shown in average ± standard error of the mean. Statistical analysis for qPCR data was performed using the doublesided Mann-Whitney U test with IBM 1 SPSS 1 Statistics version 21.0 software. A p-value < 0.05 was considered statistically significant.

Histological analysis
Full-thickness human articular cartilage specimens were harvested from the femoral heads of the non-OA and OA group, respectively. The cartilage was fixed with 4% paraformaldehyde solution in PBS and embedded in paraffin. Samples were cut into 7-μm sections and deparaffinized in xylene and dehydrated in gradually diluted ethanol from 100% to 70%. safranin-O and alcian blue/sirius red staining were performed to evaluate cartilage morphology [19].

Comparative analysis between Japanese secondary OA and Caucasian primary OA microarrays
Microarray data of Japanese hip OA were compared with the data by Xu et al [13] obtained from North European populations of UK citizens. In the study, total RNA was extracted from femoral heads of 10 #NOF (non-OA) and 9 hip OA patients in a similar manner. Gene expression profiling was assessed by Illumina Whole-genome Expression Array Human HT-12 V3 (Illumina Inc., Saffron Walden, UK). Differentially expressed genes between non-OA and OA (FC 0.5 and ! 2, p < 0.01) were extracted from Xu's study to match microarray data for comparison. The upregulated and downregulated genes were compared with our data. Furthermore, a pathway analysis using KEGG on the differentially expressed genes was conducted to assess potentially active pathways in OA chondrocytes between the Caucasian and Japanese populations.

Gene expression profiles of non-OA and OA
Microarray analysis identified 888 upregulated (FC ! 2) and 732 downregulated (FC 0.5) genes in OA chondrocytes compared to non-OA chondrocytes (p < 0.01), respectively. Roundrobin comparison of differentially expressed gene expression between 8 non-OA and 12 OA samples (96 combinations in total) revealed 352 upregulated and 159 downregulated genes in OA chondrocytes in > 80% of 96 combinations (S1 Table). Among these, expression levels of 65 genes and 12 genes were > 10 times more (FC ! 10) and less (FC 0.1) in OA compared to non-OA chondrocytes, respectively. The most differentially expressed genes between non-OA and OA are shown in Table 2

Functional enrichment analysis
The dominant function of differentially expressed genes within the 3 GO categories (cellular component, biological process, and molecular function) was evaluated. The differentially expressed genes identified 23 upregulated and 3 downregulated GO terms, respectively ( Table 3). The core enrichment genes in each significant function are listed in S2 Table. Biological process, cellular components, and molecular functions identified 6, 10, and 6 upregulated functions and 1, 0, and 2 downregulated functions, respectively. Within all categories, extracellular matrix (ECM) belonging to "cellular components" contained 61 differentially upregulated genes and was significant (FDR-q = 0.069).

Pathway analysis
Pathway analysis revealed 51 upregulated pathways in Japanese hip OA (p < 0.05) ( Table 4 and S3 Table). Overall, 47% of the genes composing each of the ECM-receptor interaction, focal adhesion, and protein digestion and absorption pathways were the differentially expressed in OA chondrocytes.  (Fig 2A-2L).

Histological analysis
The histological images of representative samples of each group are shown in Fig 3. In non-OA cartilage, a thick ECM with a smooth surface surrounds sparsely located, flat-shaped chondrocytes. In OA cartilage, round-shaped chondrocytes form clusters in eroded ECM with an irregular surface. Depletion of glycosaminoglycans and collagen fibers in the ECM of OA cartilage was demonstrated by safranin-O and alcian blue/sirius red staining, respectively.

Discussion
To the best of our knowledge, our study is the first whole-genome transcriptome analysis of articular chondrocytes in Japanese individuals with secondary hip OA. We also compared these data with transcriptome in a different ethnic origin. Few studies have reported on genome-wide transcriptome analysis of OA chondrocytes to date, presumably because extraction of high-quality nucleic acids from chondrocytes in articular cartilage is relatively difficult due to the amount of ECMs surrounding sparsely located chondrocytes [20]. Therefore, it is valuable to observe phenotypes of chondrocytes as they are in the articular cartilage. Gene expression profiling of chondrocytes in secondary hip OA yielded 888 upregulated and 732 downregulated genes in OA chondrocytes compared to non-OA in contrast to primary OA data [13] demonstrating 142 upregulated and 209 downregulated genes. In our The ratio is calculated by dividing the number of overexpressed genes in OA chondrocytes by the total number of genes contained in the pathway. Pathways overlapped in the comparative analysis between our secondary OA data and primary one [13] are shown in bold letters.
Most commonly upregulated genes in OA chondrocytes were ECM-related genes including COL2A1, COL5A2, OGN, PCOLCE, and SPARC. We found that ASPN was one of the most overexpressed genes in OA cartilage (FC = 17.05, p = 1.98E-06). Asporin, the protein encoded by ASPN, is an ECM component of the small leucine-rich proteoglycan family [21]. Asporin binds to collagen fiber to induce mineralization and inhibits TGF-β in cartilage differentiation [22]. Considering the facts above, ASPN can be potentially related to OA, particularly in patients with secondary OA. Indeed, a polymorphism in the ASPN gene is related to OA pathogenesis in the Japanese population [23]. Moreover, the genome copy number variation of ASPN can cause severe AD, a major cause of hip OA in Japan [24]. Our results support those previous studies conducted in Asian populations.  Conversely, some studies stated that ASPN does not influence OA etiology in UK Caucasians [25]. These facts indicate racial differences in genetic susceptibility of OA. As most secondary hip OA in Asians population is caused by AD [11], ASPN could be an important factor in hip OA pathology.
DPT is another gene remarkably expressed on OA chondrocytes in our study (FC = 58.28, p = 2.44E-06) that has not been reported previously as a differentially expressed gene. It is known to alter the function of TGFβ1 in endothelial cells and plays an important role in angiogenesis, fetal development, wound healing, and tumor metastasis [26]. Increased expression of DPT is observed in mesenchymal stem cells in chondrogenic differentiation [27]. Although the function of DPT in OA development and progression of OA is not clear, DPT expression is strongly inhibited by histone deacetylase (HDAC) 2 on human chondrocytes [18].
IGFBP7 was also uniquely expressed in our OA samples compared to the previous studies. This gene encodes insulin-like growth factor-binding protein (IGFBP) 7, one of the IGFBP family proteins that regulate IGF action depending on tissue-specific regulation [28]. IGFBP inhibits IGF-1, which promotes cartilage proliferation and differentiation. Increased IGFBP 3, 4, and 5 expression is observed in OA chondrocytes causing dysregulation of IGF-1 [29]. Previous studies have found the elevated expression of IGFBP-family genes in human OA cartilage than in healthy cartilage [30], however, no studies have reported of IGFBP7 expression. Our study is the first report to demonstrate elevated IGFBP7 expression in human OA chondrocytes.
During the onset and development of OA in the articular cartilage, various transcription factors (TF), e.g., HIF2A [31], RUNX2 [32], and SOX9 [33], regulate cartilage matrix-degenerating enzyme and inflammatory cytokine expression. These OA-related TFs play important roles in chondrocyte differentiation and OA progression. Meanwhile, these genes were not upregulated in OA chondrocytes in our study. Considering that expression patterns are known to vary by OA stage [34], HIF2A, RUNX2, and SOX9 genes might not have been overexpressed in end-stage OA in our study.
In the present study, 26 TF genes were overexpressed in OA versus non-OA chondrocytes. Of these, FOXO1 (forkhead-box protein O1) and KLF2 (Krüppel-like factor 2) expression was prominently elevated in OA chondrocytes. FOXO1 is widely expressed in various organs and regulates basic cellular physiological functions such as cell-cycle control, apoptosis, and glucose metabolism [35]. Akasaki et al reported FOXO1 overexpression in the middle zone of OA cartilage as well as the difference in cellular localization of FOXO1 between OA and healthy cartilage. Moreover, reduced FOXO1 expression in articular cartilage increased susceptibility of chondrocyte death under oxidative stress [36]. Further investigation is necessary to clarify the function of FOXO1 overexpression in OA chondrocytes.
KLF2, a gene encoding Krüppel-like factor 2, was the most prominently expressed TF in OA chondrocytes in our study. KLF2 is induced by fluid shear stress constituting the atherosclerosis pathway (http://genecodis.cnb.csic.es/). Central transcriptional regulator is also involved in matrix metalloproteinases (MMPs) regulation via NFκB in endothelial cells [37]. MMPs are well known catabolic genes that degrade ECM of articular cartilage [31]. Considering these facts, KLF2 can be a potential regulator of expression of MMPs in chondrocytes, and related to OA onset and/or progression. KLF2 expression was reportedly downregulated in human OA chondrocytes contrary to our results, wherein KLF2 upregulation suppressed MMP13 expression and ameliorated type II collagen degradation [38]. Although no information was provided on OA staging in Yuan's report, this discrepancy could be a result of the difference in OA grade, as expression patterns of specific genes in chondrocytes depend on the stage of OA [39]. GO analysis found that ECM components were the most significantly enriched GO term (including 61 genes, p < 0.001, FDR = 0.069). Some microarray analyses also indicated increased expression of ECM-related genes in OA cartilage [16,40]. Increased catabolism of the ECM in articular cartilage is playing a key role in the development and the progression of OA [41]. The high expression of ECM-related genes in OA chondrocytes is considered to be an "attempt" to remodel injured cartilage in response to altered cellular environment [42]. In this regard, our results support the findings in the previous studies.
Pathway analysis demonstrated the enriched ECM-receptor interaction, focal adhesion, and protein digestion and absorption pathways in the secondary OA chondrocytes. Cui et al [43] performed a meta-analysis based on 3 microarray databases for OA and healthy chondrocyte gene expression and confirmed that differentially expressed genes were enriched in ECMreceptor interaction and focal adhesion pathways. Based on enrichment map analysis for OA, they concluded that genes involved in these pathways can be essential for OA occurrence. Furthermore, these pathways affect non-traumatic necrosis of the femoral head [44] and are assumed to contribute to joint destruction. Protein digestion and absorption pathway is known to be associated with pancreatic neuroendocrine tumors and breast cancer as it is an essential protein degradation process for human nutrition homeostasis [45]. The relevance of this pathway to OA has not been clarified. In total 6 of 7 of upregulated pathways of Caucasian, primary OA were common with our Japanese, secondary data (indicated in bold in Table 4 and S3 Table), including the 3 pathways shown above. The interesting fact is that OA-related pathways are mostly common between secondary and primary hip OA, although the gene expression patterns in each pathway are unique to each other.
This study has several limitations. First, the gene expression patterns were obtained solely from the microarray analysis; therefore, the cause-and-effect relationship of each differentially expressed gene cannot be assumed. Second, the number of samples in the comparative study of primary OA was not large enough to generalize or determine the difference in expression patterns in each study. In addition, the difference of human races between primary and secondary OA in this study, can be a confounding factor. Further studies are needed to clarify the detailed network of the differentially expressed genes in hip OA chondrocytes as well as the difference between ethnic origins or etiologies of hip OA.
To summarize, to the best of our knowledge, this is the first study to investigate genomewide changes in gene expression patterns in secondary hip OA chondrocytes on the secondary hip OA in Japanese population. DPT, IGFBP7 and KLF2 were newly found genes that were prominently overexpressed in secondary hip OA chondrocytes. ECM components were the most significantly enriched GO term according to GO analysis. Pathway analysis revealed the ECM-receptor interaction pathway as an OA-related pathway that is compatible with the previous report in primary OA. Although the OA-related pathways were similar between the secondary and primary OA, the contents of the differentially expressed genes were mostly heterogeneous.
Supporting information S1