Abstract
Aims
Osteoarthritis (OA) is a common degenerative joint disease. The osteocyte transcriptome is highly relevant to osteocyte biology. This study aimed to explore the osteocyte transcriptome in subchondral bone affected by OA.
Methods
Gene expression profiles of OA subchondral bone were used to identify disease-relevant genes and signalling pathways. RNA-sequencing data of a bone loading model were used to identify the loading-responsive gene set. Weighted gene co-expression network analysis (WGCNA) was employed to develop the osteocyte mechanics-responsive gene signature.
Results
A group of 77 persistent genes that are highly relevant to extracellular matrix (ECM) biology and bone remodelling signalling were identified in OA subchondral lesions. A loading responsive gene set, including 446 principal genes, was highly enriched in OA medial tibial plateaus compared to lateral tibial plateaus. Of this gene set, a total of 223 genes were identified as the main contributors that were strongly associated with osteocyte functions and signalling pathways, such as ECM modelling, axon guidance, Hippo, Wnt, and transforming growth factor beta (TGF-β) signalling pathways. We limited the loading-responsive genes obtained via the osteocyte transcriptome signature to identify a subgroup of genes that are highly relevant to osteocytes, as the mechanics-responsive osteocyte signature in OA. Based on WGCNA, we found that this signature was highly co-expressed and identified three clusters, including early, late, and persistently responsive genes.
Conclusion
In this study, we identified the mechanics-responsive osteocyte signature in OA-lesioned subchondral bone.
Cite this article: Bone Joint Res 2022;11(6):362–370.
Article focus
-
This study explored the osteocyte transcriptome in OA subchondral bone and identified mechanics-responsive genes that contribute to subchondral bone remodelling.
Key messages
-
A group of 77 persistent genes were identified as being relevant to extracellular matrix (ECM) biology and bone remodelling signalling in OA subchondral lesions.
-
The loading responsive gene set was enriched in the OA medial tibial plateaus compared to the lateral tibial plateaus, associated with osteocyte functions and signalling pathways.
-
Three clusters, highly co-expressed in the mechanics-responsive osteocyte signature, were identified as early, late, and persistently responsive genes. The early responsive genes WNT5A and SEMA3D may be promising targets for early OA.
Strengths and limitations
-
This study elucidated mechanisms underlying osteocyte biology in subchondral bone remodelling and identified a mechanics-responsive osteocyte signature that may shed light on the blockage of mechanotransduction and help to modify early OA.
-
These findings advance the understanding of the mechanisms underlying gene-based regulation of subchondral bone mechanics and the signalling associated with it. However, the precise role of osteocytes remains elusive because the analysis was based on the gene expression profiles of bulk subchondral bone tissue rather than the isolated osteocytes.
-
This study focused exclusively on the connections between OA, subchondral bone, and the osteocyte transcriptome. The loss/gain of function assays would help to speculate on the cause-effect relationships between them. Data analysis based on multiple species may affect the reliability of this study.
Introduction
Osteoarthritis (OA) is characterized by progressive cartilage degradation, osteophyte formation, synovial inflammatory response, angiogenesis, nerve innervation, and subchondral bone sclerosis.1-4 Loading is considered a critical risk factor for OA.5 Mechanical loading exerts anti-catabolic and anabolic effects on cartilage and results in abnormal subchondral bone remodelling.6,7 Thus, reducing loading and/or blocking mechanotransduction in responsive joint tissues may help to modify OA.
Lesioned subchondral bone is an important individual feature of OA pathology. Stiffened and less pliable subchondral bone could transmit increased loads to overlying cartilage, leading to secondary cartilage damage and degeneration.8 Enhanced subchondral bone remodelling could be a vital mechanism of cartilage degradation because it happens in early OA preceding cartilage deterioration.9 Subchondral bone remodelling is characterized by bone resorption, sclerosis, and microarchitecture alterations.9 The heterogeneous outcomes of bone remodelling may be attributed to variance in spatial mechanics. Although disruption of osteoclast (Oc) and osteoblast (Ob) functions in subchondral bone lesions is highlighted,10 the role of osteocytes has been disregarded. Mechanical loading is a critical factor in OA, and osteocytes, which are mechano-sensitive11 and regulate Oc and Ob functions,12-14 are speculated to be involved in OA. Sclerostin, a key Wnt signalling inhibitor secreted by mature osteocytes, is responsive to mechanical stimulation in OA subchondral bone.15,16 Substantiating this, osteocytes showing notable morphological and phenotypic changes, such as rounded cell bodies with reduced dendrites and dysregulated expression of osteocyte markers, have been identified in OA samples.17 Although osteocyte-intrinsic transforming growth factor beta (TGF-β) signalling in perilacunar/canalicular remodelling (PLR) has been reported in OA development, the role of osteocytes in OA subchondral bone remains elusive.18 The osteocyte transcriptome has been defined as a group of genes highly relevant to osteocyte biology.19
This study investigated the osteocyte transcriptome via published gene expression profiles of OA subchondral bone lesions, to answer the following questions: 1) does the osteocyte transcriptome change in OA; 2) is the osteocyte transcriptome associated with OA subtypes; and 3) is there a group of mechanics-responsive genes in the osteocyte transcriptome that notably contribute to subchondral bone remodelling in OA?
Methods
Data processing and differential expression analysis
Gene expression profiles of GEO dataset GSE51588 were obtained from the National Centre for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO).20 The microarray dataset is based on the Agilent-026652 Whole Human Genome Microarray 4 × 44K v2 (Agilent, USA) and contains human OA (n = 20) and non-OA (n = 5) knee lateral and medial tibial plateaus (LT and MT, respectively). Differential expression analysis was performed using the R package limma (v3.11; R Foundation for Statistical Computing, Austria).21 Differentially expressed genes (DEGs) were identified using the criteria of false discovery rate (FDR) < 0.05.
Functional annotations of DEGs and gene set
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) of gene functional annotations were performed using the R package ClusterProfiler (v3.14.3; R Foundation for Statistical Computing).22 Based on hypergeometric distribution, enrichGO and enrichKEGG functions were employed to perform enrichment tests for GO terms and KEGG pathways. Gene set enrichment analysis (GSEA) was employed to determine whether a priori defined set of genes showed statistically significant differences between the groups or samples using GSEA software v4.1.0 (Broad Institute, USA).23 The estimated enrichment scores of specific gene sets of individual samples were estimated using the single-sample GSEA (ssGSEA) algorithm.24
Definition of bone loading-responsive gene set
The gene expression profiles of loading (n = 4) and non-loading (n = 4) bones were obtained from dataset GSE133212,25 which is based on the Illumina NextSeq 500 platform (Illumina, USA). The RNA-sequencing reads were aligned to the mm10 mouse genome (STAR v2.4.2a, Alexander Dobin, USA). Gene counts were derived from the number of uniquely aligned unambiguous reads (featureCount, v1.4.6, Wei Shi Lab, USA) and annotated using GENCODE (Release M27). All gene-level transcript counts were then normalized for library size and the DEGs (FDR < 0.05) were identified by comparing loading versus non-loading samples using the R package DESeq2 (v1.32.0; R Foundation for Statistical Computing),26 which was defined as the bone loading-responsive gene set.
Identification of mechanics-responsive osteocyte signature
The osteocyte transcriptome signature is defined as a group of genes (1,033 genes for humans and 1,239 genes for mice) that distinguish osteocytes from other cells.19 The candidate genes were identified by overlapping the DEGs in the comparison of MT versus LT, osteocyte transcriptome, and bone loading-responsive gene set. To identify co-expressed candidate genes, the osteocyte transcriptome signature was analyzed by weighted gene co-expression network analysis (WGCNA) using the R package WGCNA (v1.69; R Foundation for Statistical Computing)27 as previously described.28
Independent dataset validation
The independent GEO dataset GSE30322 was employed to validate the significant genes. GSE30322 contains 30 gene expression profiles of subchondral bone tissue from experimental rat OA (n = 15) and sham animals (n = 15).29
Statistical analysis
Comparisons between two groups were performed using unpaired two-tailed t-test and one-way analysis of variance (ANOVA) for comparison among more than two groups. Correlations were estimated by Pearson correlation coefficient. Pathway analysis p-values were adjusted by Benjamini & Hochberg method. All data are displayed as mean (SD) and analyzed using GraphPad Prism software version 7.0 (GraphPad, USA). Statistical significance was set at p < 0.05.
Results
Persistent genes and signalling pathways in OA subchondral bone
We summarized the clinical features of all individuals involved in this study (Supplementary Figure a). OA phenotypes assessed by the Osteoarthritis Research Society International scoring assessment system (OARSI)30 and subchondral bone parameters were only staistically significant in the comparisons of OA MT versus healthy MT and OA MT versus paired OA LT, suggesting that the MT region was more sensitive to OA. As mechanical stress is primarily transferred via MT,31 these characteristics indicated an important clue regarding the mechanics involved in OA-lesioned subchondral bone. Despite morphological parameters not being changed in OA LT samples, a distinct difference existed between OA LT and healthy LT in the human dataset GSE51588, suggesting that transcriptome shifts preceded morphological changes of subchondral bone (Supplementary Figure b). Moreover, two separated clusters were identified between OA MT and paired LT but not between healthy MT and healthy LT (Supplementary Figure b). Thus, we surmised that OA LT and OA MT samples represented the early and late OA stages, respectively, and that differences between the transcriptome of these two groups may be attributed to mechanical loading. To identify persistent genes in OA, we identified the DEGs between OA/healthy and OA MT/paired LT and overlapped these three DEG groups (Figures 1a and 1b). On comparing OA LT versus healthy LT, we identified 6,101 DEGs using FDR < 0.05. Similarly, 5,962 and 8,379 DEGs were identified in OA MT versus healthy MT and OA MT versus OA LT, respectively (Supplementary Figure c). Limited to DEGs with a robust fold-change value (log2FC < -1 or > 1), 77 common DEGs, including 52 upregulated and 25 downregulated genes, were identified as persistent genes via overlapping (Figure 1b). We employed the robust rank aggregation (RRA) method32 to identify the commonly prioritized gene list in these three comparisons. The top 25 upregulated and top ten downregulated genes are summarized (Figure 1c). Among the common genes, 19, including Aspn, Col16a1, Col2a1, Col3a1, Col5a2, Dkk3, Entpd3, Fap, Kazald1, Nid2, Ogn, Ptgfr, Ptprz1, Rhbdl2, Sema3d, Smoc2, Smpd3, Lin7a, and S100a8, were validated via an independent rat OA subchondral bone dataset (GSE30322) (Figure 1d, Supplementary Figure d, Supplementary Table i). The roles of these persistent genes in OA were determined via GO and KEGG analyses (Figure 1e). These were highly relevant to the extracellular matrix (ECM), including matrix composition, resorption, and signal transduction. ‘Axon guidance’, a pathway which is evidently relevant to osteocyte biology, was also highlighted.19
Fig. 1
Persistent genes and signalling pathways in osteoarthritis (OA) subchondral bone. a) The volcano plots of top 30 upregulated or downregulated differentially expressed genes (DEGs) in the comparisons of OA lateral tibial plateau (LT) versus healthy LT, OA medial tibial plateau (MT) versus healthy MT, and OA MT versus OA LT, respectively. b) The common genes of three comparisons, including 52 common upregulated and 25 downregulated genes. The 77 genes were defined as the persistent OA-relevant genes independent of regions. c) The heatmap plot of the persistent genes (the top 25 upregulated and top ten downregulated genes). Rank score was estimated by robust rank aggregation (RRA) method. Smaller rank score indicates greater significance in the comparison. d) Independent dataset validation of the top ten persistent genes based on the estimated rank score. OA versus sham; *p < 0.05; ** p < 0.01; all p-values estimated by unpaired two-tailed t-test. e) Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of the persistent genes. GeneRatio indicates the gene number ration in each GO term or KEGG pathway. The colour and size of each dot represents a p-value adjusted by Benjamini & Hochberg method and gene number assigned to the corresponding GO term and KEGG pathway, respectively.
Bone loading-responsive gene set is highly relevant to OA
To define the bone loading-responsive gene set, a dataset containing the gene expression profiles of four loading and four paired non-loading tibiae was recruited. Comparison between loading and non-loading groups helped to identify 2,632 genes as the transcriptome responsive to mechanics, among which 446 genes harboured robust fold changes (log2FC < -1 or > 1) (Supplementary Figure ea). These genes were highly relevant to ECM organization, focal adherence, migration, Hippo, PI3K-Akt, and Hedgehog signalling pathways (Supplementary Figure eb), suggesting a potential association between these and OA. GSEA results indicated that this gene set was significantly associated with overloading or assumed late disease stage (Figure 2a). In addition, these 446 genes were more sensitive in distinguishing between OA and healthy samples (Figure 2b). The score estimated by the loading responsive gene set (446 genes) was elevated in OA MT compared to paired OA LT, but it did not change in healthy MT compared to paired LT (Figure 2b), suggesting that the enhanced response of these genes may depend on mechanics as well as other disease-related conditions, such as inflammation.
Fig. 2
Bone loading-responsive gene set is highly relevant to osteoarthritis (OA). a) Gene Set Enrichment Analysis (GSEA) of loading responsive gene set in the osteoarthritis (OA) medial tibial plateau (MT) and the paired lateral tibial plateau (LT) samples. NES, normalized enrichment score; FDR, false discovery rate, p-values estimated by one-way analysis of variance (ANOVA) and adjusted by Benjamini-Hochberg method. Leading-edge subset was considered as the main contributing genes to the difference between the OA MT and the paired OA LT samples. b) Estimated score calculated by single-sample GSEA (ssGSEA) algorithm in the OA MT and the paired LT samples. ns, not significant; *p < 0.05; **p < 0.01; all p-values estimated by one-way analysis of variance (ANOVA). c) Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of the persistent genes. GeneRatio indicates the gene number ration in each GO term or KEGG pathway. The colour and size of each dot represents a p-value adjusted by Benjamini & Hochberg method and gene number assigned to the corresponding GO term and KEGG pathway, respectively. d) The top 50 genes in leading-edge subset (i) and the essential genes relevant to extracellular matrix (ECM) structural constituent (ii), Wnt and transforming growth factor (TGF) signalling pathways (iii), and axon guidance (iv). e) Independent dataset validation of some essential genes relevant to bone remodelling. OA versus sham; *p < 0.05; **p < 0.01; all p-values estimated by unpaired two-tailed t-test.
Among these genes, the leading-edge subset was the main contributor to the difference between OA MT and paired LT (Figure 2a). It is not surprising to find that this subset was mainly relevant to ECM and some bone modelling pathways, such as Wnt, TGF-β, Ob differentiation, and ossification (Figure 2c). We plotted the top 50 significant leading genes and the genes involved in ECM structural constitutions, axon guidance, and TGF and Wnt signalling pathways (Figure 2d). The expression of all these genes was increased in OA MT compared to paired LT. A total of 48 genes, including the top ten significant genes, COL3A1, COL18A1, CTHRC1, COL6A3, POSTN, VCAN, COL6A1, LUM, COL5A2, and ASPN, were associated with ECM structural constituents (Figure 2d). Overall 18 critical genes relevant to Wnt and TGF signalling pathways, such as TGFB3, BMP7, TGFB2, WNT5A, WNT1, WNT10B, and SOST, were identified. Among these, BMP7, WNT5A, and CAMK2B were also associated with axon guidance (Figure 2d). DCN was identified as an ECM constituent involved in the Wnt signalling pathway (Figure 2d). As these genes were relevant to OA subchondral bone remodelling, they were partially validated using an independent rat OA subchondral bone dataset (Figure 2e). After homologous mapping, we found that 60 genes were expressed in the rat validation dataset, among which 31 genes were validated and consistent with human data (Supplementary Table ii). For axon guidance, only one gene, EPHA4, was validated; for Wnt and TGF signalling, FZD1, TGFB2, TGFB3, DCN, and ACVR1 were validated; the remained 25 genes were ECM constituents, including COL3A1, COL5A2, COL6A1, COL18A1, CTHRC1, POSTN, and LUM (Figure 2f and Supplementary Table ii).
Osteocyte transcriptome signature reveals three clusters of mechanics-responsive genes
The osteocyte transcriptome signature was enriched in the OA samples compared to healthy samples (Figure 3a). The estimated enrichment scores of OA MT or LT were higher compared to those of healthy samples, while no change was observed between healthy MT and LT (Figure 3a). Cluster analysis indicated that the osteocyte signature efficiently distinguished between the two classes of OA samples (Figure 3b). This result was similar to the estimated score of the bone loading-responsive gene set (Figure 2b), suggesting a potential overlap between the osteocyte signature and the loading-responsive gene set. Indeed, we found 282 common genes by overlapping the DEGs of OA MT with paired LT, the loading-responsive gene set, and the osteocyte signature, including 277 upregulated and five downregulated genes (Figure 3c). The common genes of these three gene sets indicated an intrinsic relationship among the mechanics, osteocyte involvement, and subchondral bone lesion. The WGCNA method identified six gene modules that were groups of co-expressed and functionally relevant genes (Figure 3d). Module 5 was significantly related to disease (OA), region (MT), and the estimated scores of these two gene sets, indicating that these genes played a critical role in this module. Over 85% of the 282 common genes were assigned to this module (Supplementary Figure f). Thus, 241 co-expressed common genes were identified as osteocyte mechanics-responsive genes.
Fig. 3
Osteocyte transcriptome signature reveals three clusters of mechanics-responsive genes. a) Gene Set Enrichment Analysis (GSEA) and single-sample GSEA (ssGSEA) estimated score of loading responsive gene set. NES, normalized enrichment score; FDR, false discovery rate, p-values estimated by one-way analysis of variance (ANOVA) and adjusted by Benjamini-Hochberg method. ns, not significant; **p < 0.01. b) Principal component analysis (PCA) of the osteoarthritis (OA) medial tibial plateau (MT) and the paired lateral tibial plateau (LT) samples based on the osteocyte gene signature (Ocy). c) The overlapping of the differentially expressed genes (DEGs) (OA MT versus the paired LT), bone loading responsive genes, and Ocy. d) Heatmap of the correlation between module eigengenes and disease, region, and the estimated scores of Ocy and loading responsive gene set (LR). The value in each square reflects the Pearson correlation coefficient and the p-value in parentheses. e) The three clusters of mechanics-responsive genes. Cluster 1, the late-responsive genes; cluster 2, the early-responsive genes; cluster 3, the persistently responsive genes. The top 1 to 25 genes are labelled with red colour; the top 26 to 50 genes are labelled with orange colour. Estimated score was calculated by robust rank aggregation (RRA) algorithm to integrate the gene-disease and gene-region correlations. f) The expression of the top four genes in each cluster in the healthy and OA MT or LT samples. ns, not significant; *p < 0.05; **p < 0.01; all p-values estimated by one-way ANOVA. Blue, cluster 1; red, cluster 2; green, cluster 3. FC, fold-change value.
A gene which is significantly correlated with OA, or a related region, indicates that the revealed gene contributes to OA/healthy or MT/LT more significantly than others. As the MT region is more sensitive to OA, a gene showing good correlation with this region may display a good correlation with OA. To identify genes significantly associated with both OA and the region, we plotted these genes according to correlation between gene-disease and gene-region. However, the correlation between disease and region was still poor (r = -0.41, Pearson correlation coefficient), even though two outliers (CP and PLEKHA6) were removed from the analysis (Supplementary Figure g). This suggested that some genes regulated either the disease or the region. RRA algorithm was used to calculate a score that integrated gene-disease and gene-region correlations, with r = 0.5 as cut-off, thereby separating the genes into three clusters (Supplementary Figure h). Cluster 1, which included the genes showing a good correlation with region (r > 0.5) and a poor correlation with disease (r < 0.5), was considered to be the main contributor to the difference between the MT and LT regions (Figure 3e). This cluster remained unchanged in OA LT but significantly increased in OA MT (Figure 3f). Therefore, we defined this cluster as late-responsive genes, of which POSTN, TUBB3, GPX7, and PYCR1 were the top four, according to the integrated score of correlations. Although OA-relevant morphology phenotypes were absent in OA LT, cluster 2 genes were increased in OA LT compared to healthy LT and are thus defined as early responsive genes. Compared to cluster 2, cluster 3 was characterized by a more robust change between OA MT and paired LT and is defined as a group of persistently responsive genes. The top four genes in cluster 2 and cluster 3 included SLIT3, DPYSL3, PTPRD, and UNC5B and COL3A1, COL6A1, LUM, and GLT8D2, respectively (Figure 3f). Although WNT5A and SEMA3D in cluster 2 did not have a good rank score, they were identified as persistent.
Discussion
Osteocytes, which are terminally differentiated Obs embedded in the bone matrix, function as mechanosensors and maintain bone homeostasis and integrity.33,34 This is also applicable to subchondral bone and OA. Due to the native connection between OA and mechanics, osteocyte functions are disrupted during OA development. The role of osteocytes is important because alteration of subchondral bone during early OA may precede cartilage degeneration.35 However, the role of osteocytes in OA subchondral bone remains elusive. Investigating the role of osteocytes in OA is challenging due to difficulties in isolating osteocytes from subchondral bone caused by native embedding in the matrix and the limited number of cells available. Although immunohistochemistry and immunofluorescence may help to analyze altered gene expression, these are not convincing in some cases. More recently, success has been achieved using mouse conditional knockout models. Mazur et al36 reported that osteocytic matrix metalloproteinase 13 (MMP13) deficiency exerts a protective effect on early OA, via its effects on PLR. Bailey et al37 verified that subchondral bone remodelling in response to injury requires osteocytic TGF-β signalling. Thus, our pioneering efforts, which obtained evidence by applying accepted concepts related to bone research and OA investigations, may reveal osteocyte-related mechanisms. Although some genes and signalling pathways active in OA subchondral bone have been revealed via microarray platforms, it is challenging to apply methods such as RNA sequencing to gain an understanding of the osteocyte transcriptome.20,29 We reviewed the data of Chou et al20 and applied the osteocyte transcriptome signature defined by Youlten et al19 to identify key signalling pathways and genes relevant to osteocyte function. We found that the osteocyte transcriptome signature was associated with OA, among which the persistent and mechanics-responsive genes indicated an intrinsic relationship among the mechanics, osteocyte involvement, and subchondral bone lesion. This group of mechanics-responsive genes in the osteocyte transcriptome may contribute to subchondral bone remodelling.
Evidently, cartilage ECM is actively remodelled by chondrocytes during OA progression.38 However, shifts in the composition, arrangement, and mechanics of the ECM of OA subchondral bone are not clearly demonstrated, even though these were evidenced by micromorphological changes in subchondral bone.39 The tissue-specific composition of the osteochondral junctions of OA patients has been verified using Raman spectroscopic analysis.40 Links between biochemical shifts of ECM and cell biology, particularly osteocyte function, are rarely found. We identified OA-relevant persistent genes, suggesting dramatic ECM composition remodelling in OA subchondral bone. Collagens (COL1A2, COL2A1, COL3A1, COL5A2, COL16A1) and other ECM components, such as Asporin (ASPN) and osteonidogen (NID2), were persistently increased in OA, while S100 calcium binding protein A8 (S100A8) was decreased. S100A8 initiated early OA cartilage degradation by upregulating MMPs and aggrecanases in a mouse pre-clinical model, but its reduced expression in late OA suggested that S100A8 may not play a role in late OA.41S100A8 was decreased even in OA LTs, which showed no significant morphological progression compared with that of healthy controls. Subchondral bone is expected to respond earlier than cartilage,9 suggesting that subchondral bone S100A8 expression may act as a biomarker of early OA stage. However, the role of S100A8 in osteocytes remains unclear except for its high expression in osteocytes.41,42
The finding that a majority of the mechanical stress-responsive genes regulate ECM structural constituents indicates the link between mechanical loading and ECM remodelling. The role of some mechanical stress-responsive ECM constituents, such as periostin (POSTN), has been well demonstrated. Previous studies have shown that POSTN, mainly expressed in the periosteum and osteocytes, is increased in OA cartilage, synovial fluid, and subchondral bone.29,43,44 As a critical mechanical stress-responsive factor of secretomes, an increase in periostin upon mechanical stimulation was followed by decreased sclerostin, suggesting that POSTN may be involved in subchondral bone remodelling via sclerostin-Wnt signalling.45 We found that POSTN was a late-responsive gene, the expression of which was significantly increased only in the load-bearing MT region. Intra-articular silencing of POSTN via nanoparticle-based small interfering RNA (siRNA) treatment ameliorated subchondral bone sclerosis and heterotopic ossification,46 indicating the possibility of POSTN-targeted therapy of OA. Duan et al46 demonstrated that siRNA ameliorated the effect of OA in the MT region in mice at 18 weeks, indicating that POSTN may be a target of overloading mechanisms in late OA. It would have been helpful to know if silencing POSTN earlier in OA, such as one month after surgery, is also effective.
Since effective treatment for advanced OA is limited (except for knee arthroplasty), preventing OA initiation or progression at an early stage is vital. Animal studies have shown that subchondral bone is responsive during early OA, characterized by bone loss and microstructure changes that occur earlier than cartilage degradation.9,47,48 Thus, early responsive genes in subchondral bone, such as SLIT3, DPYSL3, PTPRD, UNC5B, WNT5A, and SEMA3D, may be promising targets. WNT5A was identified as a persistent gene that is also early responsive to mechanical stress. Our data suggest the involvement of non-canonical Wnt signalling in the mechanisms underlying mechanical loading (Figure 2d). Ob-lineage cell-derived Wnt5a is crucial for bone formation, resorption, and bone turnover.49,50 The link between OA and Wnt5a was indicated by Wnt5a being increased in OA Obs compared to normal Obs. Inhibition of Wnt5a partially attenuated abnormal mineralization and alkaline phosphatase (ALP) activity of OA Obs.51 Semaphorin 3D, encoded by SEMA3D, is a member of the semaphorin family of axon guidance factors involved in nervous system development.52 Evidently, semaphorins, such as SEMA3A, may be a potential target in OA, although the detailed roles of this family remain to be elucidated.
In conclusion, this study identified a highly co-expressing mechanics-responsive osteocyte signature in OA-lesioned subchondral bone, from which three clusters, including early, late, and persistently responsive genes, were further identified. WNT5A and SEMA3D may be potential early targets in OA subchondral bone remodelling.
References
1. Glyn-Jones S , Palmer AJR , Agricola R , et al. Osteoarthritis . Lancet . 2015 ; 386 ( 9991 ): 376 – 387 . Crossref PubMed Google Scholar
2. Chen D , Shen J , Zhao W , et al. Osteoarthritis: toward a comprehensive understanding of pathological mechanism . Bone Res . 2017 ; 5 : 16044 . Crossref PubMed Google Scholar
3. Qi X , Yu F , Wen Y , et al. Integration of transcriptome-wide association study and messenger RNA expression profile to identify genes associated with osteoarthritis . Bone Joint Res . 2020 ; 9 ( 3 ): 130 – 138 . Crossref PubMed Google Scholar
4. Duan M , Wang Q , Liu Y , Xie J . The role of TGF-β2 in cartilage development and diseases . Bone Joint Res . 2021 ; 10 ( 8 ): 474 – 487 . Crossref PubMed Google Scholar
5. Zhen G , Wen C , Jia X , et al. Inhibition of TGF-β signaling in mesenchymal stem cells of subchondral bone attenuates osteoarthritis . Nat Med . 2013 ; 19 ( 6 ): 704 – 712 . Crossref PubMed Google Scholar
6. Sun HB . Mechanical loading, cartilage degradation, and arthritis . Ann N Y Acad Sci . 2010 ; 1211 : 37 – 50 . Crossref PubMed Google Scholar
7. Chang SH , Mori D , Kobayashi H , et al. Excessive mechanical loading promotes osteoarthritis through the gremlin-1-NF-κB pathway . Nat Commun . 2019 ; 10 ( 1 ): 1442 . Crossref PubMed Google Scholar
8. Goldring MB , Goldring SR . Articular cartilage and subchondral bone in the pathogenesis of osteoarthritis . Ann N Y Acad Sci . 2010 ; 1192 : 230 – 237 . Crossref PubMed Google Scholar
9. Zhu X , Chan YT , Yung PSH , Tuan RS , Jiang Y . Subchondral bone remodeling: a therapeutic target for osteoarthritis . Front Cell Dev Biol . 2021 ; 8 : 607764 . Crossref PubMed Google Scholar
10. Maestro-Paramio L , García-Rey E , Bensiamar F , Saldaña L . Osteoblast function in patients with idiopathic osteonecrosis of the femoral head: implications for a possible novel therapy . Bone Joint Res . 2021 ; 10 ( 9 ): 619 – 628 . Crossref Google Scholar
11. Uda Y , Azab E , Sun N , Shi C , Pajevic PD . Osteocyte mechanobiology . Curr Osteoporos Rep . 2017 ; 15 ( 4 ): 318 – 325 . Crossref PubMed Google Scholar
12. Capulli M , Paone R , Rucci N . Osteoblast and osteocyte: games without frontiers . Arch Biochem Biophys . 2014 ; 561 : 3 – 12 . Crossref PubMed Google Scholar
13. Kitaura H , Marahleh A , Ohori F , et al. Osteocyte-related cytokines regulate osteoclast formation and bone resorption . Int J Mol Sci . 2020 ; 21 ( 14 ): 5169 . Crossref PubMed Google Scholar
14. Stewart S , Darwood A , Masouros S , Higgins C , Ramasamy A . Mechanotransduction in osteogenesis . Bone Joint Res . 2020 ; 9 ( 1 ): 1 – 14 . Crossref PubMed Google Scholar
15. Bouaziz W , Funck-Brentano T , Lin H , et al. Loss of sclerostin promotes osteoarthritis in mice via β-catenin-dependent and -independent Wnt pathways . Arthritis Res Ther . 2015 ; 17 ( 1 ): 24 . Crossref PubMed Google Scholar
16. Chang JC , Christiansen BA , Murugesh DK , et al. SOST/sclerostin improves posttraumatic osteoarthritis and inhibits MMP2/3 expression after injury . J Bone Miner Res . 2018 ; 33 ( 6 ): 1105 – 1113 . Crossref PubMed Google Scholar
17. Jaiprakash A , Prasadam I , Feng JQ , Liu Y , Crawford R , Xiao Y . Phenotypic characterization of osteoarthritic osteocytes from the sclerotic zones: a possible pathological role in subchondral bone sclerosis . Int J Biol Sci . 2012 ; 8 ( 3 ): 406 – 417 . Crossref PubMed Google Scholar
18. Dole NS , Mazur CM , Acevedo C , et al. Osteocyte-intrinsic TGF-β signaling regulates bone quality through perilacunar/canalicular remodeling . Cell Rep . 2017 ; 21 ( 9 ): 2585 – 2596 . Crossref PubMed Google Scholar
19. Youlten SE , Kemp JP , Logan JG , et al. Osteocyte transcriptome mapping identifies a molecular landscape controlling skeletal homeostasis and susceptibility to skeletal disease . Nat Commun . 2021 ; 12 ( 1 ): 2444 . Crossref PubMed Google Scholar
20. Chou CH , Wu CC , Song IW , et al. Genome-wide expression profiles of subchondral bone in osteoarthritis . Arthritis Res Ther . 2013 ; 15 ( 6 ): R190 . Crossref PubMed Google Scholar
21. Ritchie ME , Phipson B , Wu D , et al. limma powers differential expression analyses for RNA-sequencing and microarray studies . Nucleic Acids Res . 2015 ; 43 ( 7 ): e47 . Crossref PubMed Google Scholar
22. Yu G , Wang LG , Han Y , He QY . clusterProfiler: an R package for comparing biological themes among gene clusters . OMICS . 2012 ; 16 ( 5 ): 284 – 287 . Crossref PubMed Google Scholar
23. Subramanian A , Tamayo P , Mootha VK , et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles . Proc Natl Acad Sci U S A . 2005 ; 102 ( 43 ): 15545 – 15550 . Crossref PubMed Google Scholar
24. Hänzelmann S , Castelo R , Guinney J . GSVA: gene set variation analysis for microarray and RNA-seq data . BMC Bioinformatics . 2013 ; 14 : 7 . Crossref PubMed Google Scholar
25. Davis JL , Cox L , Shao C , et al. Conditional activation of NF-κB inducing kinase (NIK) in the osteolineage enhances both basal and loading-induced bone formation . J Bone Miner Res . 2019 ; 34 ( 11 ): 2087 – 2100 . Crossref PubMed Google Scholar
26. Love MI , Huber W , Anders S . Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 . Genome Biol . 2014 ; 15 ( 12 ): 550 . Crossref PubMed Google Scholar
27. Langfelder P , Horvath S . WGCNA: an R package for weighted correlation network analysis . BMC Bioinformatics . 2008 ; 9 : 559 . Crossref PubMed Google Scholar
28. Cheng C , Zhou J , Chen R , et al. Predicted disease-specific immune infiltration patterns decode the potential mechanisms of long non-coding RNAs in Primary Sjogren’s Syndrome . Front Immunol . 2021 ; 12 : 624614 . Crossref Google Scholar
29. Zhang R , Fang H , Chen Y , et al. Gene expression analyses of subchondral bone in early experimental osteoarthritis by microarray . PLoS One . 2012 ; 7 ( 2 ): e32356 . Crossref PubMed Google Scholar
30. Kraus VB , Huebner JL , DeGroot J , Bendele A . The OARSI histopathology initiative - recommendations for histological assessments of osteoarthritis in the guinea pig . Osteoarthr Cartil . 2010 ; 18 Suppl 3 ( Suppl 3 ): S35 - 52 . Crossref PubMed Google Scholar
31. Arjmand H , Nazemi M , Kontulainen SA , et al. Mechanical metrics of the proximal tibia are precise and differentiate osteoarthritic and normal knees: a finite element study . Sci Rep . 2018 ; 8 ( 1 ): 11478 . Crossref PubMed Google Scholar
32. Kolde R , Laur S , Adler P , Vilo J . Robust rank aggregation for gene list integration and meta-analysis . Bioinformatics . 2012 ; 28 ( 4 ): 573 – 580 . Crossref PubMed Google Scholar
33. Qin L , Liu W , Cao H , Xiao G . Molecular mechanosensors in osteocytes . Bone Res . 2020 ; 8 : 23 . Crossref PubMed Google Scholar
34. Creecy A , Damrath JG , Wallace JM . Control of bone matrix properties by osteocytes . Front Endocrinol (Lausanne) . 2021 ; 11 : 578477 . Crossref PubMed Google Scholar
35. Chen Y , Hu Y , Yu YE , et al. Subchondral trabecular rod loss and plate thickening in the development of osteoarthritis . J Bone Miner Res . 2018 ; 33 ( 2 ): 316 – 327 . Crossref PubMed Google Scholar
36. Mazur CM , Woo JJ , Yee CS , et al. Osteocyte dysfunction promotes osteoarthritis through MMP13-dependent suppression of subchondral bone homeostasis . Bone Res . 2019 ; 7 : 34 . Crossref PubMed Google Scholar
37. Bailey KN , Nguyen J , Yee CS , Dole NS , Dang A , Alliston T . Mechanosensitive control of articular cartilage and subchondral bone homeostasis in mice requires osteocytic transforming growth factor β signaling . Arthritis Rheumatol . 2021 ; 73 ( 3 ): 414 – 425 . Crossref PubMed Google Scholar
38. Maldonado M , Nam J . The role of changes in extracellular matrix of cartilage in the presence of inflammation on the pathology of osteoarthritis . Biomed Res Int . 2013 ; 2013 : 284873 . Crossref PubMed Google Scholar
39. Aho OM , Finnilä M , Thevenot J , Saarakkala S , Lehenkari P . Subchondral bone histology and grading in osteoarthritis . PLoS One . 2017 ; 12 ( 3 ): e0173726 . Crossref PubMed Google Scholar
40. Das Gupta S , Finnilä MAJ , Karhula SS , et al. Raman microspectroscopic analysis of the tissue-specific composition of the human osteochondral junction in osteoarthritis: a pilot study . Acta Biomater . 2020 ; 106 : 145 – 155 . Crossref PubMed Google Scholar
41. Zreiqat H , Belluoccio D , Smith MM , et al. S100A8 and S100A9 in experimental osteoarthritis . Arthritis Res Ther . 2010 ; 12 ( 1 ): R16 . Crossref PubMed Google Scholar
42. Zimmerman SM , Dimori M , Heard-Lipsmeyer ME , Morello R . The osteocyte transcriptome is extensively dysregulated in mouse models of osteogenesis imperfecta . JBMR Plus . 2019 ; 3 ( 7 ): e10171 . Crossref PubMed Google Scholar
43. Chijimatsu R , Kunugiza Y , Taniyama Y , Nakamura N , Tomita T , Yoshikawa H . Expression and pathological effects of periostin in human osteoarthritis cartilage . BMC Musculoskelet Disord . 2015 ; 16 : 215 . Crossref PubMed Google Scholar
44. Tajika Y , Moue T , Ishikawa S , et al. Influence of periostin on synoviocytes in knee osteoarthritis . In Vivo . 2017 ; 31 ( 1 ): 69 – 77 . Crossref PubMed Google Scholar
45. Lv J , Sun X , Ma J , et al. Involvement of periostin-sclerostin-Wnt/β-catenin signaling pathway in the prevention of neurectomy-induced bone loss by naringin . Biochem Biophys Res Commun . 2015 ; 468 ( 4 ): 587 – 593 . Crossref PubMed Google Scholar
46. Duan X , Cai L , Pham CTN , et al. Amelioration of posttraumatic osteoarthritis in mice using intraarticular silencing of periostin via nanoparticle-based small interfering RNA . Arthritis Rheumatol . 2021 ; 73 ( 12 ): 2249 – 2260 . Crossref PubMed Google Scholar
47. Hao X , Wang S , Zhang J , Xu T . Effects of body weight-supported treadmill training on cartilage-subchondral bone unit in the rat model of posttraumatic osteoarthritis . J Orthop Res . 2021 ; 39 ( 6 ): 1227 – 1235 . Crossref PubMed Google Scholar
48. Hao X , Zhang J , Shang X , et al. Exercise modifies the disease-relevant gut microbial shifts in post-traumatic osteoarthritis rats . Bone Joint Res . 2022 ; 11 ( 4 ): 214 – 225 . Crossref PubMed Google Scholar
49. Maeda K , Takahashi N , Kobayashi Y . Roles of Wnt signals in bone resorption during physiological and pathological states . J Mol Med (Berl) . 2013 ; 91 ( 1 ): 15 – 23 . Crossref PubMed Google Scholar
50. Maeda K , Kobayashi Y , Koide M , et al. The regulation of bone metabolism and disorders by wnt signaling . Int J Mol Sci . 2019 ; 20 ( 22 ): 5525 . Crossref PubMed Google Scholar
51. Martineau X , Abed É , Martel-Pelletier J , Pelletier J-P , Lajeunesse D . Alteration of Wnt5a expression and of the non-canonical Wnt/PCP and Wnt/PKC-Ca2+ pathways in human osteoarthritis osteoblasts . PLoS One . 2017 ; 12 ( 8 ): e0180711 . Crossref PubMed Google Scholar
52. Garcia S . Role of semaphorins in immunopathologies and rheumatic diseases . Int J Mol Sci . 2019 ; 20 ( 2 ): 374 . Crossref PubMed Google Scholar
Author contributions
J. Zhou: Investigation, Resources, Data curation, Writing – original draft.
Z. He: Investigation.
J. Cui: Writing – review & editing.
X. Liao: Investigation.
H. Cao: Investigation.
Y. Shibata: Investigation.
T. Miyazaki: Investigation.
J. Zhang: Conceptualization, Investigation, Methodology, Writing – original draft, Writing – review & editing, Visualization.
Funding statement
The authors disclose receipt of the following financial or material support for the research, authorship, and/or publication of this article: financial support from Grants-in-Aid for Research Activity Start-up (19K24154) from Japan Society for the Promotion of Science.
ICMJE COI statement
The authors have declared no conflict of interests regarding the publication of this paper.
Acknowledgements
The authors thank Dr Shuang Liang (Tongji hospital, Wuhan, China) for helping with the manuscript revision in this study.
Open access funding
Open access funding was provided by Grants-in-Aid for Research Activity Start-up (19K24154) from Japan Society for the Promotion of Science.
Supplementary material
Figures showing clinical features, principal component analysis, differentially expressed genes, validation, bone loading-responsive gene set, Venn plot, and correlations of the genes and disease or region. Also included are tables showing validation of persistent genes and critical bone loading-responsive genes by an independent rat osteoarthritis subchondral bone dataset GSE30322.
© 2022 Author(s) et al.This is an open-access article distributed under the terms of the Creative Commons Attribution Non-Commercial No Derivatives (CC BY-NC-ND 4.0) licence, which permits the copying and redistribution of the work only, and provided the original author and source are credited. See https://creativecommons.org/licenses/by-nc-nd/4.0/
