- Research article
- Open Access
Identification of key regulators responsible for dysregulated networks in osteoarthritis by large-scale expression analysis
Journal of Orthopaedic Surgery and Research volume 16, Article number: 259 (2021)
Osteoarthritis (OA) is a worldwide musculoskeletal disorder. However, disease-modifying therapies for OA are not available. Here, we aimed to characterize the molecular signatures of OA and to identify novel therapeutic targets and strategies to improve the treatment of OA.
We collected genome-wide transcriptome data performed on 132 OA and 74 normal human cartilage or synovium tissues from 7 independent datasets. Differential gene expression analysis and functional enrichment were performed to identify genes and pathways that were dysregulated in OA. The computational drug repurposing method was used to uncover drugs that could be repurposed to treat OA.
We identified several pathways associated with the development of OA, such as extracellular matrix organization, inflammation, bone development, and ossification. By protein-protein interaction (PPI) network analysis, we prioritized several hub genes, such as JUN, CDKN1A, VEGFA, and FOXO3. Moreover, we repurposed several FDA-approved drugs, such as cardiac glycosides, that could be used in the treatment of OA.
We proposed that the hub genes we identified would play a role in cartilage homeostasis and could be important diagnostic and therapeutic targets. Drugs such as cardiac glycosides provided new possibilities for the treatment of OA.
Osteoarthritis (OA) is the most prevalent musculoskeletal disease and a leading cause of disability in the elderly [1, 2]. OA is characterized by cartilage degeneration, the formation of osteophytes, subchondral bone remodeling, and pathological changes of the meniscus and synovitis . It causes pain, joint stiffness, and disability and leads to severe economic and social burden . Pain is the dominant symptom and a significant driver of clinical decision-making . Treatment of OA includes alleviating pain, controlling inflammation, and slowing down tissue degradation . However, there is currently no effective pharmaceutical treatment for OA that can decelerate the progression of the disease since the precise mechanisms of the pathogenesis of OA remains largely undetermined .
The integration of genome-scale transcriptomic profiling of different patient cohorts improves the understanding of molecular changes during OA progression and provides a scientific rationale for the development of novel treatment strategies. RNA-seq and microarray technology are widely used high-throughput genotyping methods that measure the expression of genes on a genome-wide scale with high accuracy and reliability [8,9,10]. Previously, it has been reported that during the pathogenesis of OA, the molecular signature of articular cartilage and synovial membrane underwent huge changes . Expression of genes involved in the molecular matrix components, cell–matrix adhesion, and ossification had dramatic changes, which typically marked different stages of endochondral ossification or transient cartilage differentiation [11,12,13]. These changes result in remodeling of the matrix and lead to the impaired function of the tissue . Notably, local low-grade inflammation of articular cartilage and synovial from patients with OA was also observed. During the progression of OA, pro-inflammatory cytokines such as IL-1, IL-6, TNF or interferon , were secreted from diseased tissues, which further damaged surrounding tissues and led to cartilage degradation . Therapeutic strategies aimed at countering inflammation is promising but does not arrest the progressive degeneration of articular cartilage [14, 16, 17]. Thus, identifying the molecular basis of OA progression and the development of novel therapeutic strategies to improve the treatment of OA remain the key issues in OA research.
To characterize OA at the molecular level and to uncover the pathogenesis mechanisms, we collected and integrated transcriptome data of hundreds of OA and normal tissues from 7 independent studies. By comparing the gene expression pattern of OA samples with normal samples, we identified significantly differentially expressed genes (DEGs) and investigated the enriched signaling pathways of these DEGs. Using protein–protein interaction (PPI) network analysis, we discovered several hub genes that function as critical regulators in OA development and might be ideal drug targets. Furthermore, we used a computational drug repurposing method to uncover several drugs that could be repurposed to treat OA.
Literature search and data collection
We used the keywords “Osteoarthritis”, “RNA-seq or microarray”, and “Dataset” on PubMed (https://www.ncbi.nlm.nih.gov/pubmed/), GEO (http://www.ncbi.nlm.nih.gov/geo/), or ArrayExpress (https://www.ebi.ac.uk/arrayexpress/) to find relevant publications or transcriptome datasets of OA. Transcriptome data of synovial or cartilage from OA patients were included in the analysis. In the data filtering and quality control process, we selected datasets with more than three replications for each condition and measuring over 10,000 genes to ensure the coverage of genes. Altogether, we collected 132 OA samples and 74 normal samples from 7 independent datasets for further analysis.
Data normalization and removal of batch effects
For the microarray data, raw “.CEL” file was downloaded. R program (version 3.5.1; http://www.r-project.org) was used for data analysis. “oligo”  and “affy”  R packages were applied for the data processing, and then the Robust Multi-array Average (RMA) method  was used for background correction, normalization, and probe set summarization. For RNA-Seq data, the sequence data were aligned to the human reference genome hg38 using STAR (v2.5.3a). RSEM (v1.2.28) was applied to map aligned reads and to generate a gene count matrix by default parameters. The expression matrix was normalized using the quantile normalization method. Residual technical batch effects arising due to heterogeneous data platforms were corrected using the ComBat  function.
Filtering of differentially expressed genes
For microarray data, “limma”  was used to perform the differential gene expression analysis of each dataset. For RNA-seq data, differential expression analysis was performed using Deseq2 . False discovery rate (FDR) was applied to carry out the correction of multiple testing using the Benjamini and Hochberg (BH) method. In this study, genes with |log2fold change (FC)| > 1.2 and FDR < 0.1 were selected as the threshold for differentially expressed genes (DEGs). Genes that were differentially expressed in over half of the datasets with the sample condition were selected and used for further analysis. Gene set enrichment analysis for DEGs was performed using “MAGeCKFlute” R package .
Integration of the PPI network and hub gene identification
The Search Tool for the Retrieval of Interacting Genes (STRING)  is a biological database for predicting protein interactions. The interactions between DEGs were evaluated using STRING, and gene sets with a combined score > 0.9 were defined as key DEGs. Subsequently, Cytoscape  (version 3.6.1; http://cytoscape.org/) was used to visualize the PPI network of the key DEGs that were identified. cytoHubba , a Cytoscape plugin, was used to extract the hub genes. The central elements were ranked by betweenness.
Identification of putative target genes of JUN
We searched public ChIP-seq data of JUN on the Cistrome Data Browser (http://cistrome.org/db) , a website that collected and integrated thousands of public ChIP-seq data of both human and mouse. We identified and downloaded 25 processed ChIP-seq data of JUN on Cistrome Data Browser. Genes with high RP scores had high likelihoods to be the target genes of JUN. We selected genes with mean RP scores > 1 as candidate target genes of JUN. For further filtering of these candidate genes, we performed co-expression analysis to identify genes that were co-expressed with JUN. Genes with absolute correlation value > 0.5 in each dataset were selected as co-expressed genes. Candidate target genes that were co-expressed in over half of the expression datasets were defined as putative target genes of JUN.
Drug repurposing analysis
Genes that were differentially expressed between OA and normal from both synovial and cartilage tissues were used for drug repurposing analysis. The probe name of the upregulated and downregulated genes were used as input on the cMap (https://portals.broadinstitute.org/cmap/) website. The “quick query” mode of the cMap algorithm was used for the drug repurposing analysis.
Statistical analyses were performed using the R software (http://www.R-project.org/). Statistical analyses gathering more than two groups were performed using ANOVA. Otherwise, for two groups, statistical analyses were performed using the unpaired t test. Multiple hypothesis testing corrections were applied where multiple hypotheses were tested and were indicated by the use of FDR.
Outline of data analysis
The objective of this study was to use differentially expressed genes (DEGs) to identify dysregulated genes and pathways of OA and to uncover potential novel pharmacological strategies. To achieve this purpose, we collected hundreds of transcriptome data of patients with OA and normal tissues from 7 independent datasets (Table 1). The subsequent analyses focused on prioritizing key regulators by protein–protein interaction (PPI) network and identifying potent drugs that can be repurposed to treat OA by the computational drug repurposing method (Fig. 1).
Molecular characteristic of OA
To get a list of OA-related DEGs, we compared gene expression profiles of OA patients with those of normal tissues. We collected transcriptome datasets from synovial or cartilage tissues of OA patients and healthy donors. To distinguish DEGs between OA and normal, we performed the differential analysis of each dataset. Among the four transcriptome datasets of synovial tissues, there were 579 DEGs from GSE1919, 2441 DEGs from GSE12021, 2431 DEGs from GSE55235, and 2423 DEGs from GSE55457. The overlap of the DEGs between each dataset was shown in Fig. 2a. Genes that were differentially expressed in two of the four datasets were selected for gene set enrichment analysis. Genes associated with the lysosome, oxidative phosphorylation, extracellular matrix organization, endopeptidase activity, skeletal system development, and collagen-containing extracellular matrix were upregulated in synovial tissues of OA patients than healthy donors (Fig. 2b). In contrast, compared with normal tissues, genes related to IL-17 signaling pathway, circadian rhythm, positive regulation of p38mapk cascade, NOD-like receptor signaling pathway, Foxo signaling pathway, cellular senescence, negative regulation of erk1 and erk2 cascade pathway, and cellular response to hypoxia were downregulated in patients with OA (Fig. 2b).
To distinguish whether the cartilage tissues had different gene expression patterns from what we have observed in synovial tissues, we also collected three transcriptome datasets from cartilage tissues. By comparing the gene expression pattern between patients with OA and healthy donors, we obtained 6307 DEGs from E-MTAB-6266, 6037 DEGs from GSE114007, and 1578 DEGs from GSE117999. The overlap of the DEGs among these three datasets was shown in Fig. 2c. Genes that were differentially expressed in two of the three datasets were selected for further analysis. Gene ontology analysis of the DEGs revealed that genes involved in collagen-containing extracellular matrix, collagen binding, N-acetyglucosamine metabolic process, collagen fibril organization, skeletal system development, ossification, and osteoblast differentiation were upregulated in OA samples (Fig. 2d); while biological processes, including mTOR signaling pathway, Foxo signaling pathway, regulation of cell proliferation, cellular senescence, positive regulation of p38mapk cascade, and ERBB signaling pathway, were downregulated in OA samples (Fig. 2d). Although there were uniquely enriched pathways for synovial and cartilage tissues, we also found some pathways that were enriched in both synovial and cartilage tissues. Such pathways might play essential roles in the development and pathology of OA.
To better understand the differences between OA and normal on the molecular expression pattern, we focused on the genes that were differentially expressed in both synovial and cartilage tissues for further analysis. Among the 3036 DEGs of cartilage and 2620 DEGs of synovial, there were 555 genes differentially expressed in both tissues (Fig. 3a, Supplementary table 1). Next, we performed functional enrichment analysis for these 555 common DEGs. Compared with normal tissues, genes involved in glycine, serine and threonine metabolism, extracellular matrix organization, regulation of T cell activation, peptidase activator activity, bone development, ossification, and collagen-containing extracellular matrix were enriched in OA samples (Fig. 3b), while pathways related to negative regulation of erk1 and erk2 cascade, cellular senescence, response to hypoxia, circadian rhythm, positive regulation of p38mapk cascade, and HIF-1 signaling pathway were downregulated in OA samples (Fig. 3c). Genes involved in these pathways might be responsible for the cartilage homeostasis and the pathology of OA.
Protein–protein interaction network of DEGs
To discover potential critical regulators of OA from DEGs, we performed protein–protein interaction (PPI) network analysis (Fig. 4a). We used cytoHubba, a Cytoscape plugin, to extract the hub genes in the network. The top 10 hub genes were JUN, ITGB1, VEGFA, CDKN1A, CDK2, PSMD3, CUL2, FOXO3, SGK1, and INSR (Fig.4b). These hub genes might associate with the development and pathogenesis of OA. JUN is a transcription factor that plays an essential role in regulating cellular proliferation and apoptosis. JUN is in the center of the network and interacts with a large number of DEGs in the network, suggesting the crucial role of JUN in regulating the expression of these DEGs. Integrins are major surface receptors of chondrocytes, and integrin β1 (ITGB1) has been proven to suppress chondrocyte hypertrophy and accelerate chondrocyte proliferation by in vitro studies and mouse models [34,35,36,37], suggesting that ITGB1 is involved in OA progression and OA-induced cartilage degradation. VEGF-A is the founding member of the VEGF family and is classically referred to as VEGF . As reported, the expression levels of VEGF were associated with OA progression and OA-specific pathologies, such as cartilage degeneration, osteophyte formation, synovitis, and pain. Moreover, a wide range of studies suggested that inhibition of VEGF signaling reduces OA progression and benefits patients with OA . CDKN1A encodes a potent cyclin-dependent kinase inhibitor and functions as a regulator of cell cycle progression at G1. Shinsuke et al. reported that CDKN1A-deficient mice were more susceptible to OA-related changes in vivo, suggesting that CDKN1A modulation might constitute a possible therapeutic strategy for OA treatment . FOXO3 belongs to the forkhead family of transcription factors, which are characterized by a distinct forkhead domain . FOXO3 functions as a trigger for apoptosis by regulating the expression of genes necessary for cell death . Akasaki et al. reported that normal articular cartilage had a tissue-specific signature of FOXO1 and FOXO3 but not FOXO4 proteins. In OA cartilage, chondrocytes showed altered FOXO activation, which suggested FOXO might play a role in OA progression . SGK1 belongs to the serine/threonine-protein kinase subfamily, which contributes to the regulation of a wide variety of physiological activities, such as membrane transports, cell growth, proliferation, and apoptosis . Knockdown of SGK1 alleviates the IL-1β-induced chondrocyte anabolic and catabolic imbalance by activating autophagy in human chondrocytes . Downregulation of SGK1 attenuates OA progression . The biological role of these hub genes in the development and progression of OA had been studied previously, suggesting that targeting these hub genes might be an optimal novel treatment for OA. However, more experiments were needed to better elucidate the mechanism of action of these hub genes.
JUN functions as a key TF that associated with the development of OA
Since JUN was in the central of the network (Fig. 4a), we focused on JUN for further analysis. We found that most genes interacting with JUN directly in the PPI network were differentially expressed between OA and normal (Fig. 5a and b), suggesting the importance of JUN in regulating the expression of these DEGs. To identify the potential target genes that were regulated by JUN, we downloaded the processed ChIP-seq data of JUN on the Cistrome Data Browser (http://cistrome.org/db) . Cistrome Data Browser collected thousands of human and mouse samples with ChIP-seq data and facilitated searches for putative target genes of transcription factors by regulatory potential model . Regulatory potential (RP) score for each gene was calculated, reflecting the likelihood of the transcription factor being a direct regulator of that gene, and genes with high RP scores were putative targets of the transcription factor. By analyzing 25 public ChIP-seq data of JUN, we selected 3250 genes (Supplementary table 2) with mean RP scores > 1 as candidate target genes of JUN. For further filtering of the candidate genes, we performed co-expression analysis. Among the 3250 candidate targets, 214 genes had high correlations with JUN, which were selected for further analysis (Supplementary table 3). Functional enrichment analysis of the selected 214 genes revealed that pathways including the ATF-2 transcription factor network (part of the AP-1 complex), AP-1 transcription factor network, TGF-β signaling pathway, osteoclast differentiation, and ERBB1 downstream signaling were enriched (Fig. 5c). JUN was in the AP-1 protein family, and the top 2 enriched pathways were associated with the AP-1 family network, which suggested the creditability of our data analysis method. Among these, the TGF-β signaling pathway has been reported to play a critical role in the development and progression of OA by driving chondrocytes toward hypertrophy, promoting osteoprogenitor cell differentiation into osteoblasts, mediating synthesis of cartilage-specific extracellular matrix components, and angiogenesis in subchondral bone . The effects of TGF-β on the modulation of extracellular matrix components were dependent on the activation of JNK (c-Jun N-terminal Kinase), which in turn modulates the activity of c-Jun . Thus, JUN might serve as a critical regulator of the development of OA by regulating the activity of TGF-β.
Drug repurposing for OA
Since there are no effective interventions to decelerate the progression of OA now, drug repurposing of the Food and Drug Administration (FDA)-approved therapeutic agents is a particularly attractive approach to improve the treatment of OA. Computational techniques for predictive repurposing offer a relatively efficient and authentic method of identifying testable hypotheses that may be translated into the clinic . The Connectivity Map (cMap), which was established by the Broad Institute, consists of gene expression data generated by dosing of more than 1300 compounds in hundreds of cell lines . cMap (https://portals.broadinstitute.org/cmap/) has been successfully used to make drug repurposing predictions for a number of disease conditions . We applied cMap to explore potential drug repositioning opportunities for OA. Among the top 15 listed drugs, MG-262 ranked first (Table 2). MG-262 is a potent proteasome inhibitor that selectively and reversibly inhibits the chymotryptic activity of the proteasome [52, 53]. The proteasome inhibitors have shown anti-inflammatory activities in the animal models of arthritis, psoriasis, colitis, and other inflammatory conditions . Some studies discovered that proteasome inhibitors promoted bone growth in a cell-based screen . Although the biological function of MG-262 in the treatment of OA still remains an unexploited field, this analysis provided new possibilities for the treatment of OA. Of note, among the top 15 drugs, there were four drugs classified as cardiac glycosides, such as ouabain and digoxin, suggesting the potential clinical implications of cardiac glycosides in OA. Cardiac glycosides function by inhibiting the Na+/K+ ATPase (NKA) . In addition, cardiac glycosides have been found to decrease inflammatory symptoms in different animal models of acute and chronic inflammation . Thus, cardiac glycosides might have therapeutic benefit in the treatment of OA by countering the inflammation-induced articular cartilage degradation. Since the safety of FDA-approved drugs has been sufficiently verified, drug repurposing for approved drugs offers a less risky and more rapid potential for the investment of novel therapeutic strategies.
Here, we used a large-scale data integration method to characterize the molecular signatures of OA, which extended our understanding of the disease mechanisms. We also identified several essential regulators of OA, such as JUN, VEGFA, and FOXO3, which might provide a scientific rationale for the development of novel pharmacological therapies. Of note, we found JUN, a crucial dysregulated transcription factor, plays a central role in regulating the aberrant gene expression pattern in OA. Moreover, we used a computational drug repurposing method to identify potential FDA-approved drugs that can be repurposed to improve the treatment of OA.
In this study, we integrated and analyzed multiple published datasets. We were able to find genes that were consistently differentially expressed between OA and normal among different datasets and different tissues. Our analysis confirmed some findings from previous studies using genome-wide gene expression analyses. Common findings of these studies are the differential expression of genes involved in matrix-degrading enzymes (MMPs, ADAMTS), collagen organization, and inflammation [7, 58]. In our analysis, we also found increased expression of genes associated with collagen-containing extracellular matrix (DPT, GPC4, COL8A2, FBLN5, COL3A1), ossification (CDH11, COL1A1, SPARC, KAZALD1), and bone development (DYM, PTGER4, SPARC, SULF) in patients with OA, which is suggestive of active remodeling of cartilage homeostasis during OA pathogenesis. During the early stages of OA, the molecular composition and organization of the extracellular matrix are altered first . The articular chondrocytes exhibit increased cell proliferation and matrix synthesis for the purpose of initiating repairing for pathological injury [59, 60]. Changes in the composition and structure of the articular cartilage further stimulate chondrocytes to produce more catabolic factors involved in cartilage degradation. Thus, the expression of genes involved in the carbohydrate metabolism and extracellular matrix components were upregulated. Moreover, chronic low-grade inflammation has also been found to contribute to the development and progression of OA . During OA progression, the entire synovial joints were involved in the inflammation process . Pro-inflammatory factors, such as IL-1β and TNF-α, as well as chemokines, were reported to contribute to the systemic inflammation that led to the activation of NF-κB signaling in both synovial cells and chondrocytes . Based on these studies, multiple novel pharmacological strategies have emerged, including anti-inflammatory mediators (anti-IL-1 , anti-TNF-α , and anti-IL-6 ) and inhibition of catabolic pathways (Wnt, ADAMTS, and cathepsin K) . Apart from these findings, we were able to find some essential regulators that might associate with the pathogenesis of OA, such as ITGB1, PSMD3, CUL2, and INSR. Since the association of these hub genes with the development of OA has been reported, they could serve as important diagnostic and/or therapeutic targets for OA. However, more mechanistic interpretation and clinical trial data are needed to clarify the efficacy in the treatment of OA.
Since synovial and cartilage play different roles in the bone homeostasis and pathology of OA, we compared the gene expression pattern between OA and normal using transcriptome datasets of each tissue separately. This analysis identified some tissue-specific DEGs between OA and normal. Gene sets associated with N-acetylglucosamine metabolic process, glycosaminoglycan biosynthesis and metabolism, NF-KB transcription factor activity, ERBB signaling pathway, P53 signaling pathway, and mTOR signaling pathway were uniquely differentially expressed in cartilage tissues. While biological pathways involved in lysosome, oxidative phosphorylation, endopeptidase activity, Nod-like receptor signaling pathway, and IL-17 signaling pathway were uniquely differentially expressed in synovial tissues. In the development of OA, the pathological changes of synovial and cartilage were different. The histological pattern of synovium in OA patients is characterized by synovial lining hyperplasia, increased vascularity, and sublining fibrosis  with the phenotypic shift of chondrocytes, including surface fibrillation, degradation of cartilage matrix, chondrocyte clusters appearance, and vascular penetration from the subchondral bone . The different phenotypic changes of synovial and cartilage might explain the different OA-related DEGs between these two tissues.
Our study also identified JUN, a transcription factor, as a key regulator of these DEGs. JUN is one of the members of the Activator protein 1 (AP-1) family proteins . AP-1 family proteins are basic leucine zipper (bZIP) transcription factors that consist of Jun (c-Jun, JunB, and JunD), Fos (c-Fos, FosB, Fra-1, and Fra-2), Jun dimerization partners (JDP1 and JDP2), and the closely related activating transcription factors (ATF2, LRF1/ ATF3, and B-ATF) subfamilies . AP-1 family proteins are implicated in the regulation of a variety of cellular processes, including proliferation and survival, differentiation, apoptosis, cell migration, and transformation . JUN has been reported to play a crucial role in regulating cell proliferation and apoptosis . Ventura et al. reported that the c-Jun NH2-terminal kinase JNK signaling pathway contributes to the regulation of TGF-β-mediated biological responses . TGF-β is crucial for cartilage maintenance, and the lack of TGF-β results in OA-like changes . Thus, JUN might regulate the development of OA by coordinating with TGF-β signaling.
Drug discovery is a time-consuming, laborious, costly, and high-risk process. Drug repositioning is efficient, economical, and low-risk compared with the traditional drug development process. There are multiple drug repurposing methods generated, including computational approaches, biological experimental approaches, and mixed approaches. Among these, computational drug repurposing is the most powerful and the most widely used. Previously, Soul et al. constructed a data portal that provides an exploration and comparison platform for analyzed skeletal transcriptomics data . This data portal also provides data repurposing mode using L1000 data, one of the cMap project resources. Using a computational drug repurposing method, we found cardiac glycosides might be repurposed in the treatment of OA. Cardiac glycosides are drugs that inhibit the Na+/K+ ATPases and are applied to treat heart failure and certain irregular heartbeats. However, recent studies reported that cardiac glycosides are a novel class of broad-spectrum senolytics for therapeutic applications in many age-related disorders , including osteoarthritis. Cardiac glycosides were capable of reducing the number of senescent cells, diminishing the level of local inflammation, and improved some metabolic and physical fitness parameters that decline with aging in some animal models [77, 78]. Drug repurposing can be highly attractive as a potentially cheaper and faster route to market. However, successful drug repositioning requires a deep understanding of biological mechanisms—from known overlaps of mechanisms to re-innovation of a new molecule, to find a new mechanism, dosing, route of administration, and new target. Although our drug repurposing analysis offers novel and valuable options for developing strategies to treat OA, the effectiveness of the predicted drugs in treating OA, such as cardiac glycosides, needs more systemic and detailed validation.
There are several limitations of our studies. First, osteoarthritis is typically described as a heterogeneous disease with complex pathogenesis. Different patients might have different mechanistic pathways, such that the mechanisms of OA in the elderly might be different from those after a joint injury in a younger adult or in obese individuals. The molecular signatures we discovered in this study by transcriptome data analysis might not be representative in some subtypes of OA patients. Integrated analysis of multi-omics data, including epigenetics (DNA methylation, histone post-translational modification, and/or non-coding RNA), metabolomics, and proteomics (by LC-MS) data, are needed to better elucidate the heterogeneity of OA. Such detailed and systematic analyses offer important mechanistic and potentially therapeutic insights into OA. Finally, we have highlighted several essential regulators and uncovered that cardiac glycosides might benefit the patients with OA by countering the inflammation. However, further experimental studies are needed to validate their clinical implications.
Our study integrated multiple public transcriptome data sets together, which provides more comprehensive and reliable insights into the genetic alterations associated with the disease phenotype. Moreover, we used a computational drug repurposing method to identify potent drug candidates to improve the treatment of OA.
In summary, we used bioinformatics analysis to identify a group of differentially expressed genes between OA and normal tissues. By protein–protein interaction network analysis, we identified several osteoarthritis-related hub genes, which might be potential diagnostic or therapeutic targets for osteoarthritis. In addition, we found JUN, a transcription factor, functions as a crucial regulator by the analysis of public ChIP-seq data and co-expression analysis. Moreover, we used a computational drug repurposing method to identify potent drugs that can be repurposed to treat OA.
Availability of data and materials
The datasets used and/or analyzed in this study are available from the corresponding author on reasonable request.
Differentially expressed genes
Robust Multi-array Average
Search Tool for the Retrieval of Interacting Genes
Woolf AD, Pfleger B. Burden of major musculoskeletal conditions. Bull World Health Organ. 2003;81(9):646–56.
Chou CH, Wu CC, Song IW, Chuang HP, Lu LS, Chang JH, et al. Genome-wide expression profiles of subchondral bone in osteoarthritis. Arthritis Res Ther. 2013;15(6):R190. https://doi.org/10.1186/ar4380.
Katz JN. Lumbar disc disorders and low-back pain: socioeconomic factors and consequences. J Bone Joint Surg Am. 2006;88(Suppl 2):21–4.
Shepherd C, Reese AE, Reynard LN, Loughlin J. Expression analysis of the osteoarthritis genetic susceptibility mapping to the matrix Gla protein gene MGP. Arthritis Res Ther. 2019;21(1):149. https://doi.org/10.1186/s13075-019-1934-7.
Glyn-Jones S, Palmer AJ, Agricola R, Price AJ, Vincent TL, Weinans H, et al. Osteoarthritis. Lancet. 2015;386(9991):376–87. https://doi.org/10.1016/S0140-6736(14)60802-3.
Bijlsma JW, Berenbaum F, Lafeber FP. Osteoarthritis: an update with relevance for clinical practice. Lancet. 2011;377(9783):2115–26. https://doi.org/10.1016/S0140-6736(11)60243-2.
Malaise O, de Seny D. Therapeutic advances in arthritis diseases. Biochem Pharmacol. 2019;165:1–3. https://doi.org/10.1016/j.bcp.2019.04.014.
Chen G, Shi T, Shi L. Characterizing and annotating the genome using RNA-seq data. Sci China Life Sci. 2017;60(2):116–25. https://doi.org/10.1007/s11427-015-0349-4.
Nagalakshmi U, Waern K, Snyder M. RNA-Seq: a method for comprehensive transcriptome analysis. Curr Protoc Mol Biol. 2010;Chapter 4:Unit 4 11 1-13.
Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5(7):621–8. https://doi.org/10.1038/nmeth.1226.
Pitsillides AA, Beier F. Cartilage biology in osteoarthritis--lessons from developmental biology. Nat Rev Rheumatol. 2011;7(11):654–63. https://doi.org/10.1038/nrrheum.2011.129.
Xu YK, Ke Y, Wang B, Lin JH. The role of MCP-1-CCR2 ligand-receptor axis in chondrocyte degradation and disease progress in knee osteoarthritis. Biol Res. 2015;48(1):64. https://doi.org/10.1186/s40659-015-0057-0.
Heinegard D, Saxne T. The role of the cartilage matrix in osteoarthritis. Nat Rev Rheumatol. 2011;7(1):50–6. https://doi.org/10.1038/nrrheum.2010.198.
Jaswal AP, Bandyopadhyay A. Re-examining osteoarthritis therapy from a developmental biologist's perspective. Biochem Pharmacol. 2019;165:17–23. https://doi.org/10.1016/j.bcp.2019.03.020.
Henrotin Y, Lambert C, Richette P. Importance of synovitis in osteoarthritis: evidence for the use of glycosaminoglycans against synovial inflammation. Semin Arthritis Rheum. 2014;43(5):579–87. https://doi.org/10.1016/j.semarthrit.2013.10.005.
Rannou F, Poiraudeau S. Non-pharmacological approaches for the treatment of osteoarthritis. Best Pract Res Clin Rheumatol. 2010;24(1):93–106. https://doi.org/10.1016/j.berh.2009.08.013.
Zhang W, Moskowitz RW, Nuki G, Abramson S, Altman RD, Arden N, et al. OARSI recommendations for the management of hip and knee osteoarthritis, Part II: OARSI evidence-based, expert consensus guidelines. Osteoarthritis Cartilage. 2008;16(2):137–62. https://doi.org/10.1016/j.joca.2007.12.013.
Carvalho BS, Irizarry RA. A framework for oligonucleotide microarray preprocessing. Bioinformatics. 2010;26(19):2363–7. https://doi.org/10.1093/bioinformatics/btq431.
Gautier L, Cope L, Bolstad BM, Irizarry RA. affy--analysis of Affymetrix GeneChip data at the probe level. Bioinformatics. 2004;20(3):307–15.
Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, et al. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003;4(2):249–64. https://doi.org/10.1093/biostatistics/4.2.249.
Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012;28(6):882–3. https://doi.org/10.1093/bioinformatics/bts034.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.
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. https://doi.org/10.1186/s13059-014-0550-8.
Wang B, Wang M, Zhang W, Xiao T, Chen CH, Wu A, et al. Integrative analysis of pooled CRISPR genetic screens using MAGeCKFlute. Nat Protoc. 2019;14(3):756–80. https://doi.org/10.1038/s41596-018-0113-7.
Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47(D1):D607–13. https://doi.org/10.1093/nar/gky1131.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504. https://doi.org/10.1101/gr.1239303.
Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8 Suppl 4:S11.
Zheng R, Wan C, Mei S, Qin Q, Wu Q, Sun H, et al. Cistrome Data Browser: expanded datasets and new tools for gene regulatory analysis. Nucleic Acids Res. 2019;47(D1):D729–35. https://doi.org/10.1093/nar/gky1094.
Woetzel D, Huber R, Kupfer P, Pohlers D, Pfaff M, Driesch D, et al. Identification of rheumatoid arthritis and osteoarthritis patients by transcriptome-based rule set generation. Arthritis Res Ther. 2014;16(2):R84. https://doi.org/10.1186/ar4526.
Huber R, Hummert C, Gausmann U, Pohlers D, Koczan D, Guthke R, et al. Identification of intra-group, inter-individual, and gene-specific variances in mRNA expression profiles in the rheumatoid arthritis synovial membrane. Arthritis Res Ther. 2008;10(4):R98. https://doi.org/10.1186/ar2485.
Ungethuem U, Haeupl T, Witt H, Koczan D, Krenn V, Huber H, et al. Molecular signatures and new candidates to target the pathogenesis of rheumatoid arthritis. Physiol Genomics. 2010;42A(4):267–82. https://doi.org/10.1152/physiolgenomics.00004.2010.
Fisch KM, Gamini R, Alvarez-Garcia O, Akagi R, Saito M, Muramatsu Y, et al. Identification of transcription factors responsible for dysregulated networks in human osteoarthritis cartilage by global gene expression analysis. Osteoarthritis Cartilage. 2018;26(11):1531–8. https://doi.org/10.1016/j.joca.2018.07.012.
Soul J, Dunn SL, Anand S, Serracino-Inglott F, Schwartz JM, Boot-Handford RP, et al. Stratification of knee osteoarthritis: two major patient subgroups identified by genome-wide expression analysis of articular cartilage. Ann Rheum Dis. 2018;77(3):423. https://doi.org/10.1136/annrheumdis-2017-212603.
Enomoto M, Leboy PS, Menko AS, Boettiger D. Beta 1 integrins mediate chondrocyte interaction with type I collagen, type II collagen, and fibronectin. Exp Cell Res. 1993;205(2):276–85. https://doi.org/10.1006/excr.1993.1087.
Raducanu A, Hunziker EB, Drosse I, Aszodi A. Beta1 integrin deficiency results in multiple abnormalities of the knee joint. J Biol Chem. 2009;284(35):23780–92. https://doi.org/10.1074/jbc.M109.039347.
Bengtsson T, Aszodi A, Nicolae C, Hunziker EB, Lundgren-Akerlund E, Fassler R. Loss of alpha10beta1 integrin expression leads to moderate dysfunction of growth plate chondrocytes. J Cell Sci. 2005;118(Pt 5):929–36. https://doi.org/10.1242/jcs.01678.
Lian C, Wang X, Qiu X, Wu Z, Gao B, Liu L, et al. Collagen type II suppresses articular chondrocyte hypertrophy and osteoarthritis progression by promoting integrin beta1-SMAD1 interaction. Bone Res. 2019;7(1):8. https://doi.org/10.1038/s41413-019-0046-y.
Ferrara N, Gerber HP, LeCouter J. The biology of VEGF and its receptors. Nat Med. 2003;9(6):669–76. https://doi.org/10.1038/nm0603-669.
Hamilton JL, Nagao M, Levine BR, Chen D, Olsen BR, Im HJ. Targeting VEGF and its receptors for the treatment of osteoarthritis and associated pain. J Bone Miner Res. 2016;31(5):911–24. https://doi.org/10.1002/jbmr.2828.
Kihara S, Hayashi S, Hashimoto S, Kanzaki N, Takayama K, Matsumoto T, et al. Cyclin-dependent kinase inhibitor-1-deficient mice are susceptible to osteoarthritis associated with enhanced inflammation. J Bone Miner Res. 2017;32(5):991–1001. https://doi.org/10.1002/jbmr.3080.
Stefanetti RJ, Voisin S, Russell A, Lamon S. Recent advances in understanding the role of FOXO3, F1000Res 7; 2018.
Ekoff M, Kaufmann T, Engstrom M, Motoyama N, Villunger A, Jonsson JI, et al. The BH3-only protein Puma plays an essential role in cytokine deprivation induced apoptosis of mast cells. Blood. 2007;110(9):3209–17. https://doi.org/10.1182/blood-2007-02-073957.
Akasaki Y, Hasegawa A, Saito M, Asahara H, Iwamoto Y, Lotz MK. Dysregulated FOXO transcription factors in articular cartilage in aging and osteoarthritis. Osteoarthritis Cartilage. 2014;22(1):162–70. https://doi.org/10.1016/j.joca.2013.11.004.
Lang F, Strutz-Seebohm N, Seebohm G, Lang UE. Significance of SGK1 in the regulation of neuronal function. J Physiol. 2010;588(Pt 18):3349–54. https://doi.org/10.1113/jphysiol.2010.190926.
Huang W, Cheng C, Shan WS, Ding ZF, Liu FE, Lu W, et al. Knockdown of SGK1 alleviates the IL-1beta-induced chondrocyte anabolic and catabolic imbalance by activating FoxO1-mediated autophagy in human chondrocytes. FEBS J. 2020;287(1):94–107. https://doi.org/10.1111/febs.15009.
Wang Z, Ni S, Zhang H, Fan Y, Xia L, Li N. Silencing SGK1 alleviates osteoarthritis through epigenetic regulation of CREB1 and ABCA1 expression. Life Sci. 2021;268:118733. https://doi.org/10.1016/j.lfs.2020.118733.
Tang Q, Chen Y, Meyer C, Geistlinger T, Lupien M, Wang Q, et al. A comprehensive view of nuclear receptor cancer cistromes. Cancer Res. 2011;71(22):6940–7. https://doi.org/10.1158/0008-5472.CAN-11-2091.
Shen J, Li S, Chen D. TGF-beta signaling and the development of osteoarthritis, Bone Res 2; 2014.
Hocevar BA, Brown TL, Howe PH. TGF-beta induces fibronectin synthesis through a c-Jun N-terminal kinase-dependent, Smad4-independent pathway. EMBO J. 1999;18(5):1345–56. https://doi.org/10.1093/emboj/18.5.1345.
Pushpakom S, Iorio F, Eyers PA, Escott KJ, Hopper S, Wells A, et al. Drug repurposing: progress, challenges and recommendations. Nat Rev Drug Discov. 2019;18(1):41–58. https://doi.org/10.1038/nrd.2018.168.
Lamb J, Crawford ED, Peck D, Modell JW, Blat IC, Wrobel MJ, et al. The Connectivity Map: using gene-expression signatures to connect small molecules, genes, and disease. Science. 2006;313(5795):1929–35. https://doi.org/10.1126/science.1132939.
Pujols L, Fernandez-Bertolin L, Fuentes-Prado M, Alobid I, Roca-Ferrer J, Agell N, et al. Proteasome inhibition reduces proliferation, collagen expression, and inflammatory cytokine production in nasal mucosa and polyp fibroblasts. J Pharmacol Exp Ther. 2012;343(1):184–97. https://doi.org/10.1124/jpet.111.190710.
Frase H, Hudak J, Lee I. Identification of the proteasome inhibitor MG262 as a potent ATP-dependent inhibitor of the Salmonella enterica serovar Typhimurium Lon protease. Biochemistry. 2006;45(27):8264–74. https://doi.org/10.1021/bi060542e.
Kisselev AF, van der Linden WA, Overkleeft HS. Proteasome inhibitors: an expanding army attacking a unique target. Chem Biol. 2012;19(1):99–115. https://doi.org/10.1016/j.chembiol.2012.01.003.
Garrett IR, Chen D, Gutierrez G, Zhao M, Escobedo A, Rossini G, et al. Selective inhibitors of the osteoblast proteasome stimulate bone formation in vivo and in vitro. J Clin Invest. 2003;111(11):1771–82. https://doi.org/10.1172/JCI16198.
Koren G, Soldin SJ. Cardiac glycosides. Clin Lab Med. 1987;7(3):587–606. https://doi.org/10.1016/S0272-2712(18)30733-9.
Furst R, Zundorf I, Dingermann T. New knowledge about old drugs: the anti-inflammatory properties of cardiac glycosides. Planta Med. 2017;83(12-13):977–84. https://doi.org/10.1055/s-0043-105390.
Li H, Yang HH, Sun ZG, Tang HB, Min JK. Whole-transcriptome sequencing of knee joint cartilage from osteoarthritis patients. Bone Joint Res. 2019;8(7):290–303. https://doi.org/10.1302/2046-3758.87.BJR-2018-0297.R1.
Lee AS, Ellman MB, Yan D, Kroin JS, Cole BJ, van Wijnen AJ, et al. A current review of molecular mechanisms regarding osteoarthritis and pain. Gene. 2013;527(2):440–7. https://doi.org/10.1016/j.gene.2013.05.069.
Im HJ, Li X, Muddasani P, Kim GH, Davis F, Rangan J, et al. Basic fibroblast growth factor accelerates matrix degradation via a neuro-endocrine pathway in human adult articular chondrocytes. J Cell Physiol. 2008;215(2):452–63. https://doi.org/10.1002/jcp.21317.
Shen J, Abu-Amer Y, O'Keefe RJ, McAlinden A. Inflammation and epigenetic regulation in osteoarthritis. Connect Tissue Res. 2017;58(1):49–63. https://doi.org/10.1080/03008207.2016.1208655.
Goldring MB, Otero M. Inflammation in osteoarthritis. Curr Opin Rheumatol. 2011;23(5):471–8. https://doi.org/10.1097/BOR.0b013e328349c2b1.
Chen D, Shen J, Zhao W, Wang T, Han L, Hamilton JL, et al. Osteoarthritis: toward a comprehensive understanding of pathological mechanism. Bone Res. 2017;5(1):16044. https://doi.org/10.1038/boneres.2016.44.
Singer JW, Fleischman A, Al-Fayoumi S, Mascarenhas JO, Yu Q, Agarwal A. Inhibition of interleukin-1 receptor-associated kinase 1 (IRAK1) as a therapeutic strategy. Oncotarget. 2018;9(70):33416–39. https://doi.org/10.18632/oncotarget.26058.
Grunke M, Schulze-Koops H. Successful treatment of inflammatory knee osteoarthritis with tumour necrosis factor blockade. Ann Rheum Dis. 2006;65(4):555–6. https://doi.org/10.1136/ard.2006.053272.
Brenn D, Richter F, Schaible HG. Sensitization of unmyelinated sensory fibers of the joint nerve to mechanical stimuli by interleukin-6 in the rat: an inflammatory mechanism of joint pain. Arthritis Rheum. 2007;56(1):351–9. https://doi.org/10.1002/art.22282.
Pap T, Bertrand J. Syndecans in cartilage breakdown and synovial inflammation. Nat Rev Rheumatol. 2013;9(1):43–55. https://doi.org/10.1038/nrrheum.2012.178.
Mathiessen A, Conaghan PG. Synovitis in osteoarthritis: current understanding with therapeutic implications. Arthritis Res Ther. 2017;19(1):18. https://doi.org/10.1186/s13075-017-1229-9.
Goldring MB. Articular cartilage degradation in osteoarthritis. HSS J. 2012;8(1):7–9. https://doi.org/10.1007/s11420-011-9250-z.
Karin M, Liu Z, Zandi E. AP-1 function and regulation. Curr Opin Cell Biol. 1997;9(2):240–6. https://doi.org/10.1016/S0955-0674(97)80068-3.
Vesely PW, Staber PB, Hoefler G, Kenner L. Translational regulation mechanisms of AP-1 proteins. Mutat Res. 2009;682(1):7–12. https://doi.org/10.1016/j.mrrev.2009.01.001.
Shaulian E, Karin M. AP-1 in cell proliferation and survival. Oncogene. 2001;20(19):2390–400. https://doi.org/10.1038/sj.onc.1204383.
Kappelmann M, Bosserhoff A, Kuphal S. AP-1/c-Jun transcription factors: regulation and function in malignant melanoma. Eur J Cell Biol. 2014;93(1-2):76–81. https://doi.org/10.1016/j.ejcb.2013.10.003.
Ventura JJ, Kennedy NJ, Flavell RA, Davis RJ. JNK regulates autocrine expression of TGF-beta1. Mol Cell. 2004;15(2):269–78. https://doi.org/10.1016/j.molcel.2004.06.007.
Soul J, Hardingham TE, Boot-Handford RP, Schwartz JM. SkeletalVis: an exploration and meta-analysis data portal of cross-species skeletal transcriptomics data. Bioinformatics. 2019;35(13):2283–90. https://doi.org/10.1093/bioinformatics/bty947.
Martin N, Soriani O, Bernard D. Cardiac glycosides as senolytic compounds. Trends Mol Med. 2020;26(3):243–5. https://doi.org/10.1016/j.molmed.2020.01.001.
He S, Sharpless NE. Senescence in health and disease. Cell. 2017;169(6):1000–11. https://doi.org/10.1016/j.cell.2017.05.015.
Paez-Ribes M, Gonzalez-Gualda E, Doherty GJ, Munoz-Espin D. Targeting senescent cells in translational medicine. EMBO Mol Med. 2019;11(12):e10234. https://doi.org/10.15252/emmm.201810234.
The authors wish to thank the Gene Expression Omnibus (GEO) and ArrayExpress database for making their processed data publicly available.
This research was partially supported by the National Natural Science Foundation of China (82002282).
Ethics approval and consent to participate
Consent for publication
All the authors approved the manuscript for publication.
The authors declare that they have no conflicts of interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Shi, S., Wan, F., Zhou, Z. et al. Identification of key regulators responsible for dysregulated networks in osteoarthritis by large-scale expression analysis. J Orthop Surg Res 16, 259 (2021). https://doi.org/10.1186/s13018-021-02402-9
- Transcriptional profiling
- Drug repurposing