Identification of Key Regulators Responsible for Dysregulated Networks in Osteoarthritis by Large-Scale Expression Analysis

doi:10.21203/rs.3.rs-47411/v1

Download PDF

Research article

Identification of Key Regulators Responsible for Dysregulated Networks in Osteoarthritis by Large-Scale Expression Analysis

https://doi.org/10.21203/rs.3.rs-47411/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 14 Apr, 2021

Read the published version in Journal of Orthopaedic Surgery and Research →

You are reading this older preprint version

Read the latest preprint version →

Background: 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.

Methods: We collected genome-wide transcriptome data performed on 161 OA and 88 normal human cartilage or synovium tissues from 9 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.

Results: We identified several pathways associated with the development of OA, such as extracellular matrix adhesion, inflammation, carbohydrate metabolism, and collagen organization. By protein-protein interaction (PPI) network analysis, we prioritized several hub genes, such as JUN, IL-6, VAMP8, and SOCS3. Moreover, we repurposed several FDA-approved drugs, such as cardiac glycosides, that could be used in the treatment of OA.

Conclusions: We proposed that the hub genes we identified will play a role in cartilage homeostasis and could be important diagnostic and therapeutic targets. Drugs such as protein kinase inhibitors and cardiac glycosides could improve the treatment of OA.

Orthopedics

Osteoarthritis

transcriptional profiling

drug repurposing

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[3]. It causes pain, joint stiffness, disability, and leads to severe economic and social burden[4]. Pain is the dominant symptom and a significant driver of clinical decision-making[5]. Treatment of OA includes pain alleviation, controlling inflammation, and slowing down tissue degradation[6]. 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[7].

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–10]. Previously, it has been reported that during the pathogenesis of OA, the molecular signature of articular cartilage and synovial membrane underwent huge changes[11]. 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–13]. These changes result in remodeling of the matrix and lead to the impaired function of the tissue[14]. 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[15], were secreted from diseased tissues, which further damaged surrounding tissues and led to cartilage degradation[7]. 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 remains 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 9 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”, “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 167 OA samples and 88 normal samples from 9 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”[18] and “affy”[19] R packages were applied for the data processing, and then the Robust Multi-array Average (RMA) method[20] was used for background correction, normalization, and probe-set summarization. For RNA-Seq data, the sequence data were aligned to the human reference genome using STAR. RSEM was applied to map aligned reads and to generate a gene count matrix. 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[21] function.

Filtering of differentially expressed genes

For microarray data, “limma”[22] was used to perform the differential gene expression analysis of each dataset. For RNA-Seq data, differential expression analysis was performed using Deseq2[23]. 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[24].

Integration of the PPI network and hub gene identification

The Search Tool for the Retrieval of Interacting Genes (STRING)[25] 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[26] (version 3.6.1; http://cytoscape.org/) was used to visualize the PPI network of the key DEGs that were identified. cytoHubba[27], 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)[28], a website 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.

Statistical 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 9 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).

Table 1

Information about the collected datasets
	1	2	3	4	5	6	7	8	9
Data	GSE55457	GSE55235	GSE12021	GSE32317	GSE1919	GSE82107	GSE117999	GSE114007	E-MTAB-6266
Type	Array	Array	Array	Array	Array	Array	Array	RNA-Seq	RNA-Seq
Tissue	Synorial	Synorial	Synorial	Synorial	Synorial	Synorial	Cartilage	Cartilage	Cartilage
Normal	10	10	9	7	5	7	12	18	10
OA	10	10	10	19	5	10	12	20	65

Table 1

top 10 drugs predicted by cMap
rank	cmap name	mean	n	enrichment	specificity	Description
1	H-7	0.717	4	0.998	0	protein kinase inhibitor, Inhibiting the activity of PKA, PKG, PKC and MLCK
2	anisomycin	-0.843	4	-0.993	0	antibiotic, inhibiting eukaryotic protein synthesis
3	MG-262	-0.722	3	-0.983	0	proteasome inhibitor, inhibiting the chymotryptic activity of the proteasome
4	lanatoside C	-0.767	6	-0.982	0	cardiac glycosides, inhibiting the plasma membrane Na+/K + ATPase
5	ouabain	-0.767	4	-0.977	0.0088	cardiac glycoside, inhibiting the Na+/K + ATPase
6	digoxin	-0.699	4	-0.968	0	cardiac glycoside, inhibiting the Na+/K + ATPase
7	digoxigenin	-0.645	5	-0.952	0	derivative of the cardiac glycoside digoxin
8	helveticoside	-0.665	6	-0.95	0	digitalis-like cardiac glycosidic
9	phenoxybenzamine	-0.576	4	-0.943	0.0091	alpha-adrenoceptor antagonist,used as an anti-hypertensive
10	digitoxigenin	-0.696	4	-0.942	0	the aglycone of digitoxin

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. Genes that were differentially expressed in 5 of the 9 datasets were selected as DEGs for further analysis (Supplementary table 1). To identify the biological processes for these DEGs, we performed functional enrichment analysis. Compared with normal tissues, genes involved in the lysosome, collagen-containing extracellular matrix, collagen binding, ossification, cell-matrix adhesion, integrin-mediated signaling pathway, and cytokine production were up-regulated in OA samples (Fig. 2A). While biological processes, including Foxo signaling pathway, skeletal muscle cell differentiation, positive regulation of p38mapk cascade, negative regulation of ERK1 and ERK2 cascade, Hif-1 signaling pathway, regulation of cell cycle, and Erbb signaling pathway were enriched in normal tissues (Fig. 2B). 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. 3A). We used cytoHubba, a Cytoscape plugin, to extract the hub genes in the network. The top 10 hub genes were JUN, POLR2L, IL6, SOCS3, VAMP8, GAA, AGA, EGR1, ICAM, and PTAFR (Fig. 3B). These hub genes might be responsible for 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. IL-6, another hub gene, has been reported to be responsible for the pathogenesis of OA[29]. IL-6 pathway regulates the activation of JAK/STAT signaling and elevates the activity of MMP, ADAMs, and ADAMTS, leading to cartilage extracellular matrix protein (ECM) degradation and the loss of accessory matrix proteins[29]. Thus, IL-6 could be an important target for future drug development in OA. As another hub gene, VAMP8 encoded proteins for the fusion of cellular membranes and was reported to mediate lysosomal exocytosis in hypertrophic chondrocytes in the context of OA[30]. Another hub gene is SOCS3, which is a potent regulator of cytokine signaling. Activation of SOCS3 in bone cells, including osteoclasts, chondrocytes, and osteoblasts, has been reported[31]. SOCS3 could serve as a key signaling molecule in bone cell-mediated inflammatory responses[32]. These hub genes might be ideal therapeutic targets for experimental biologists to validate their role in the development of OA.

JUN functions as a key TF that associated with the development of OA

Since JUN was in the central of the network (Fig. 3A), 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. 4A), suggesting the importance of JUN in regulating the expression of these DEGs. To identify 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) [28]. 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 [33]. 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. 4B). 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[34]. 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[35]. Thus, JUN might serve as a critical regulator of the development of OA by regulating the activity of TGF‐β.

Identification and characterization of molecular subgroups

Osteoarthritis is typically described as a multifaceted and heterogeneous syndrome characterized by variable clinical features, biochemical/genetic characteristics, and responses to treatments[36]. Each of the common osteoarthritis risk factors might instigate a different mechanistic pathway leading to osteoarthritis. A number of stratifications have been proposed on the basis of specific pathological processes to classify different mechanistic subgroups, which include an increased inflammatory component, mechanical overload, metabolic alterations, and cell senescence[37]. Knowledge of the heterogeneity of OA would help better define potential treatment targets, improve the ability to select appropriate patients for future clinical interventions and adapt outcome measures to the situation[38]. To identify potential subgroups based on their expression profiles, we performed KNN clustering. OA patients were divided into two subgroups based on the expression of differentially expressed genes between OA and normal (Fig. 5A). Segregation of OA subgroups was also reproduced by principal component analysis (Fig. 5B). To identify characteristic molecular signaling pathways enriched in each subgroup, we performed functional enrichment analyses for each subgroup's differentially expressed genes. Pathways related to inflammatory response, such as neutrophil degranulation, MHC-II antigen presentation, co-stimulation by the CD28 family, and interferon-gamma signaling, were enriched in subgroup1 (Fig. 5C). Patients in subgroup2 had increased expression levels of genes associated with extracellular matrix organization, integrin cell surface interactions, and ECM proteoglycans (Fig. 5D). The different expression patterns of OA subgroups might explain the heterogeneity of response to OA treatment.

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 for improving 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[39]. The Connectivity Map (cMap), which was established by the Broad Institute, consists of gene expression data generated by dosing of more than 1,300 compounds in hundreds of cell lines[40]. cMap (https://portals.broadinstitute.org/cmap/) has been successfully used to make drug repurposing predictions for a number of disease conditions[39]. We applied cMap to explore potential drug repositioning opportunities for OA. Among the top 10 listed drugs, H-7 ranked first. H-7 is a protein kinase inhibitor, which inhibits the activity of PKA, PKG, PKC, and MLCK (Table 2). It has been reported that PKC was up-regulated in osteoarthritic cartilage, and overexpression PKC prevented cytokine-mediated NF-kB activation in primary chondrocytes[41]. PKC could also function as a necessary component of the IL-1 and TNF signaling mediated catabolic destruction[41]. Drugs inhibiting PKC, such as H-7, provides new possibilities for the treatment of OA. Of note, among the top 10 drugs, there were six drugs classified as cardiac glycosides, such as lanatoside C, ouabain, and digoxin, suggesting the potential clinical implications of cardiac glycosides in OA. Cardiac glycosides function by inhibiting the Na+/K + ATPase (NKA)[42]. In addition, cardiac glycosides have been found to decrease inflammatory symptoms in different animal models of acute and chronic inflammation[43]. Thus, Cardiac glycosides might have therapeutic benefit in the treatment of OA by countering the inflammation-induced articular cartilage degradation. Drug repurposing for approved drugs offers a less risky and more rapid potential for the investment of novel therapeutic strategies.

Here, we used 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, IL-6, and SOCS3, 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, which confirms 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, 44]. In our analysis, we also found increased expression of genes associated with carbohydrate metabolism (BPGM, CHST2, GNPDA1, GUSB, HYAL1), integrin mediated signaling (FUT8, ITGA5, ITGB5), and ossification (COL5A2, COL6A1, CTGF, EXT2) 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[45]. The articular chondrocytes exhibit increased cell proliferation and matrix synthesis for the purpose of initiating repairing for pathological injury[45, 46]. 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 up-regulated. Moreover, chronic low-grade inflammation has also been found contributing to the development and progression of OA[47]. During OA progression, the entire synovial joints were involved in the inflammation process[48]. 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[49]. Based on these studies, multiple novel pharmacological strategies have emerged, including anti-inflammatory mediators (anti-IL-1[50], anti-TNF-α[51], and anti-IL-6[52]) and inhibition of catabolic pathways (Wnt, ADAMTS, and cathepsin K)[53]. Apart from these findings, we were able to find some essential regulators that might associate with the pathogenesis of OA. The hub genes we identified, such as IL-6, VAMP8, and SOCS3, could serve as important diagnostic and/or therapeutic targets for OA.

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[54]. 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[55]. 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[56]. JUN has been reported to play a crucial role in regulating cell proliferation and apoptosis[57]. Ventura et al. reported that the c-Jun NH2-terminal kinase JNK signaling pathway contributes to the regulation of TGF-β-mediated biological responses[58]. TGF-β is crucial for cartilage maintenance, and lack of TGF-β results in OA-like changes[58]. Thus, JUN might regulate the development of OA by coordinating with TGF-β signaling.

In our analysis, we found two subgroups of OA patients. Patients in subgroup 1 had increased expression levels of genes related with inflammatory response, such as IL1B, TNF, TLR1, TLR7. Currently, countering inflammation is the mainstay of osteoarthritis therapy today. However, only a minority of patients respond to such therapy, suggesting the necessity to select appropriate patient subgroups. Patients with high expression of genes involved in inflammatory response might show a better response to anti-inflammatory drugs, such as anti-TNF and anti-IL1B. Patients in subgroup 2 had a higher expression level of extracellular matrix organization-related genes, such as MMP2, MMP28, ADAMTS5. Matrix metallopeptidases (MMPs), capable of degrading all kinds of extracellular matrix proteins, are therapeutic targets that showed efficacy in preventing joint destruction in some preclinical evaluations[59]. However, many MMPs inhibitors failed phase I or phase II clinical trials due to the limited efficacy or drug toxicity. We hypothesized that patients in subgroup 2 who had a higher expression level of extracellular matrix proteins might better respond to MMPs inhibitors. Our analysis provided new possibilities for different treatment strategies according to the gene expression pattern of different OA subgroups.

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 applied to treat heart failure and certain irregular heartbeats. However, recent studies reported that cardiac glycosides as a novel class of broad-spectrum senolytics for therapeutic applications in many age-related disorders[60], 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[61, 62]. Although the effectiveness of cardiac glycosides in treating OA needs more systemic and detailed validation, our analysis offers a novel and valuable option for patients with OA.

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. Further detailed analyses that take these co-variants (age, obesity and/or injury) into consideration may help better identify OA subgroups, which might explain the heterogeneity of OA patients and the different response to the treatment. Second, although we have discovered some molecular signatures of OA by transcriptome data analysis, multi-omics investigations, such as the integrated analysis of epigenetic (DNA methylation, histone post-translational modification, and/or noncoding RNA), transcriptomic (via RNA-seq or microarray), and proteomic (by LC-MS) changes are needed to better elucidate the changes during the pathogenesis 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 risk 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 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 osteoarthritis.

OA: Osteoarthritis; DEGs: Differentially expressed genes; PPI: Protein–protein interaction; Robust Multi-array Average (RMA); ECM: Extracellular matrix; GO: Gene Ontology; STRING: Search Tool for the Retrieval of Interacting Genes; Regulatory potential (RP).

Acknowledgements

The authors wish to thank the Gene Expression omnibus (GEO) and ArrayExpress database for making their processed data publicly available.

Author Contributions

S.S, F.W, Y.L, and F.L were involved in the study design. S.S, F.W, Z.Z, R.T, Y.L, and M.J participated in the data collection. S.S and F.W performed all the data analyses. S.S, F.W, F.L and Y.L drafted the manuscript. All authors contributed to the discussion and revision of the manuscript.

Funding

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sector.

Availability of data and materials

The datasets used and/or analyzed in this study are available from the corresponding author on reasonable request.

Ethics approval and consent to participate.

Not applicable.

Consent for publication

Not applicable.

Competing interests

All authors declare no competing interests.

A.D. Woolf, B. Pfleger, Burden of major musculoskeletal conditions, Bull World Health Organ 81(9) (2003) 646-56.
C.H. Chou, C.C. Wu, I.W. Song, H.P. Chuang, L.S. Lu, J.H. Chang, S.Y. Kuo, C.H. Lee, J.Y. Wu, Y.T. Chen, V.B. Kraus, M.T. Lee, Genome-wide expression profiles of subchondral bone in osteoarthritis, Arthritis Res Ther 15(6) (2013) R190.
J.N. Katz, Lumbar disc disorders and low-back pain: socioeconomic factors and consequences, J Bone Joint Surg Am 88 Suppl 2 (2006) 21-4.
C. Shepherd, A.E. Reese, L.N. Reynard, J. Loughlin, Expression analysis of the osteoarthritis genetic susceptibility mapping to the matrix Gla protein gene MGP, Arthritis Res Ther 21(1) (2019) 149.
S. Glyn-Jones, A.J. Palmer, R. Agricola, A.J. Price, T.L. Vincent, H. Weinans, A.J. Carr, Osteoarthritis, Lancet 386(9991) (2015) 376-87.
J.W. Bijlsma, F. Berenbaum, F.P. Lafeber, Osteoarthritis: an update with relevance for clinical practice, Lancet 377(9783) (2011) 2115-26.
O. Malaise, D. de Seny, Therapeutic advances in arthritis diseases, Biochem Pharmacol 165 (2019) 1-3.
G. Chen, T. Shi, L. Shi, Characterizing and annotating the genome using RNA-seq data, Sci China Life Sci 60(2) (2017) 116-125.
U. Nagalakshmi, K. Waern, M. Snyder, RNA-Seq: a method for comprehensive transcriptome analysis, Curr Protoc Mol Biol Chapter 4 (2010) Unit 4 11 1-13.
A. Mortazavi, B.A. Williams, K. McCue, L. Schaeffer, B. Wold, Mapping and quantifying mammalian transcriptomes by RNA-Seq, Nat Methods 5(7) (2008) 621-8.
A.A. Pitsillides, F. Beier, Cartilage biology in osteoarthritis--lessons from developmental biology, Nat Rev Rheumatol 7(11) (2011) 654-63.
Y.K. Xu, Y. Ke, B. Wang, J.H. Lin, The role of MCP-1-CCR2 ligand-receptor axis in chondrocyte degradation and disease progress in knee osteoarthritis, Biol Res 48 (2015) 64.
D. Heinegard, T. Saxne, The role of the cartilage matrix in osteoarthritis, Nat Rev Rheumatol 7(1) (2011) 50-6.
A.P. Jaswal, A. Bandyopadhyay, Re-examining osteoarthritis therapy from a developmental biologist's perspective, Biochem Pharmacol 165 (2019) 17-23.
Y. Henrotin, C. Lambert, P. Richette, Importance of synovitis in osteoarthritis: evidence for the use of glycosaminoglycans against synovial inflammation, Semin Arthritis Rheum 43(5) (2014) 579-87.
F. Rannou, S. Poiraudeau, Non-pharmacological approaches for the treatment of osteoarthritis, Best Pract Res Clin Rheumatol 24(1) (2010) 93-106.
W. Zhang, R.W. Moskowitz, G. Nuki, S. Abramson, R.D. Altman, N. Arden, S. Bierma-Zeinstra, K.D. Brandt, P. Croft, M. Doherty, M. Dougados, M. Hochberg, D.J. Hunter, K. Kwoh, L.S. Lohmander, P. Tugwell, OARSI recommendations for the management of hip and knee osteoarthritis, Part II: OARSI evidence-based, expert consensus guidelines, Osteoarthritis Cartilage 16(2) (2008) 137-62.
B.S. Carvalho, R.A. Irizarry, A framework for oligonucleotide microarray preprocessing, Bioinformatics 26(19) (2010) 2363-7.
L. Gautier, L. Cope, B.M. Bolstad, R.A. Irizarry, affy--analysis of Affymetrix GeneChip data at the probe level, Bioinformatics 20(3) (2004) 307-15.
R.A. Irizarry, B. Hobbs, F. Collin, Y.D. Beazer-Barclay, K.J. Antonellis, U. Scherf, T.P. Speed, Exploration, normalization, and summaries of high density oligonucleotide array probe level data, Biostatistics 4(2) (2003) 249-64.
J.T. Leek, W.E. Johnson, H.S. Parker, A.E. Jaffe, J.D. Storey, The sva package for removing batch effects and other unwanted variation in high-throughput experiments, Bioinformatics 28(6) (2012) 882-3.
M.E. Ritchie, B. Phipson, D. Wu, Y. Hu, C.W. Law, W. Shi, G.K. Smyth, limma powers differential expression analyses for RNA-sequencing and microarray studies, Nucleic Acids Res 43(7) (2015) e47.
M.I. Love, W. Huber, S. Anders, Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2, Genome Biol 15(12) (2014) 550.
B. Wang, M. Wang, W. Zhang, T. Xiao, C.H. Chen, A. Wu, F. Wu, N. Traugh, X. Wang, Z. Li, S. Mei, Y. Cui, S. Shi, J.J. Lipp, M. Hinterndorfer, J. Zuber, M. Brown, W. Li, X.S. Liu, Integrative analysis of pooled CRISPR genetic screens using MAGeCKFlute, Nat Protoc 14(3) (2019) 756-780.
D. Szklarczyk, A.L. Gable, D. Lyon, A. Junge, S. Wyder, J. Huerta-Cepas, M. Simonovic, N.T. Doncheva, J.H. Morris, P. Bork, L.J. Jensen, C.V. Mering, STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets, Nucleic Acids Res 47(D1) (2019) D607-D613.
P. Shannon, A. Markiel, O. Ozier, N.S. Baliga, J.T. Wang, D. Ramage, N. Amin, B. Schwikowski, T. Ideker, Cytoscape: a software environment for integrated models of biomolecular interaction networks, Genome Res 13(11) (2003) 2498-504.
C.H. Chin, S.H. Chen, H.H. Wu, C.W. Ho, M.T. Ko, C.Y. Lin, cytoHubba: identifying hub objects and sub-networks from complex interactome, BMC Syst Biol 8 Suppl 4 (2014) S11.
R. Zheng, C. Wan, S. Mei, Q. Qin, Q. Wu, H. Sun, C.H. Chen, M. Brown, X. Zhang, C.A. Meyer, X.S. Liu, Cistrome Data Browser: expanded datasets and new tools for gene regulatory analysis, Nucleic Acids Res 47(D1) (2019) D729-D735.
M. Laavola, T. Leppanen, M. Hamalainen, K. Vuolteenaho, T. Moilanen, R. Nieminen, E. Moilanen, IL-6 in Osteoarthritis: Effects of Pine Stilbenoids, Molecules 24(1) (2018).
E.R. Bastow, K. Last, S. Golub, J.L. Stow, A.C. Stanley, A.J. Fosang, Evidence for lysosomal exocytosis and release of aggrecan-degrading hydrolases from hypertrophic chondrocytes, in vitro and in vivo, Biol Open 1(4) (2012) 318-28.
T. Gui, B.S. He, Q. Gan, C. Yang, Enhanced SOCS3 in osteoarthiritis may limit both proliferation and inflammation, Biotech Histochem 92(2) (2017) 107-114.
A. Gao, T.E. Van Dyke, Role of suppressors of cytokine signaling 3 in bone inflammatory responses, Front Immunol 4 (2014) 506.
Q. Tang, Y. Chen, C. Meyer, T. Geistlinger, M. Lupien, Q. Wang, T. Liu, Y. Zhang, M. Brown, X.S. Liu, A comprehensive view of nuclear receptor cancer cistromes, Cancer Res 71(22) (2011) 6940-7.
J. Shen, S. Li, D. Chen, TGF-beta signaling and the development of osteoarthritis, Bone Res 2 (2014).
B.A. Hocevar, T.L. Brown, P.H. Howe, TGF-beta induces fibronectin synthesis through a c-Jun N-terminal kinase-dependent, Smad4-independent pathway, EMBO J 18(5) (1999) 1345-56.
L.A. Deveza, R.F. Loeser, Is osteoarthritis one disease or a collection of many?, Rheumatology (Oxford) 57(suppl_4) (2018) iv34-iv42.
S. Grassel, A. Aszodi, Osteoarthritis and Cartilage Regeneration: Focus on Pathophysiology and Molecular Mechanisms, Int J Mol Sci 20(24) (2019).
W.E. Van Spil, O. Kubassova, M. Boesen, A.C. Bay-Jensen, A. Mobasheri, Osteoarthritis phenotypes and novel therapeutic targets, Biochem Pharmacol 165 (2019) 41-48.
S. Pushpakom, F. Iorio, P.A. Eyers, K.J. Escott, S. Hopper, A. Wells, A. Doig, T. Guilliams, J. Latimer, C. McNamee, A. Norris, P. Sanseau, D. Cavalla, M. Pirmohamed, Drug repurposing: progress, challenges and recommendations, Nat Rev Drug Discov 18(1) (2019) 41-58.
J. Lamb, E.D. Crawford, D. Peck, J.W. Modell, I.C. Blat, M.J. Wrobel, J. Lerner, J.P. Brunet, A. Subramanian, K.N. Ross, M. Reich, H. Hieronymus, G. Wei, S.A. Armstrong, S.J. Haggarty, P.A. Clemons, R. Wei, S.A. Carr, E.S. Lander, T.R. Golub, The Connectivity Map: using gene-expression signatures to connect small molecules, genes, and disease, Science 313(5795) (2006) 1929-35.
E.R. LaVallie, P.S. Chockalingam, L.A. Collins-Racie, B.A. Freeman, C.C. Keohan, M. Leitges, A.J. Dorner, E.A. Morris, M.K. Majumdar, M. Arai, Protein kinase Czeta is up-regulated in osteoarthritic cartilage and is required for activation of NF-kappaB by tumor necrosis factor and interleukin-1 in articular chondrocytes, J Biol Chem 281(34) (2006) 24124-37.
G. Koren, S.J. Soldin, Cardiac glycosides, Clin Lab Med 7(3) (1987) 587-606.
R. Furst, I. Zundorf, T. Dingermann, New Knowledge About Old Drugs: The Anti-Inflammatory Properties of Cardiac Glycosides, Planta Med 83(12-13) (2017) 977-984.
H. Li, H.H. Yang, Z.G. Sun, H.B. Tang, J.K. Min, Whole-transcriptome sequencing of knee joint cartilage from osteoarthritis patients, Bone Joint Res 8(7) (2019) 290-303.
A.S. Lee, M.B. Ellman, D. Yan, J.S. Kroin, B.J. Cole, A.J. van Wijnen, H.J. Im, A current review of molecular mechanisms regarding osteoarthritis and pain, Gene 527(2) (2013) 440-7.
H.J. Im, X. Li, P. Muddasani, G.H. Kim, F. Davis, J. Rangan, C.B. Forsyth, M. Ellman, E.J. Thonar, Basic fibroblast growth factor accelerates matrix degradation via a neuro-endocrine pathway in human adult articular chondrocytes, J Cell Physiol 215(2) (2008) 452-63.
J. Shen, Y. Abu-Amer, R.J. O'Keefe, A. McAlinden, Inflammation and epigenetic regulation in osteoarthritis, Connect Tissue Res 58(1) (2017) 49-63.
M.B. Goldring, M. Otero, Inflammation in osteoarthritis, Curr Opin Rheumatol 23(5) (2011) 471-8.
D. Chen, J. Shen, W. Zhao, T. Wang, L. Han, J.L. Hamilton, H.J. Im, Osteoarthritis: toward a comprehensive understanding of pathological mechanism, Bone Res 5 (2017) 16044.
J.W. Singer, A. Fleischman, S. Al-Fayoumi, J.O. Mascarenhas, Q. Yu, A. Agarwal, Inhibition of interleukin-1 receptor-associated kinase 1 (IRAK1) as a therapeutic strategy, Oncotarget 9(70) (2018) 33416-33439.
M. Grunke, H. Schulze-Koops, Successful treatment of inflammatory knee osteoarthritis with tumour necrosis factor blockade, Ann Rheum Dis 65(4) (2006) 555-6.
D. Brenn, F. Richter, H.G. Schaible, 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 56(1) (2007) 351-9.
T. Pap, J. Bertrand, Syndecans in cartilage breakdown and synovial inflammation, Nat Rev Rheumatol 9(1) (2013) 43-55.
M. Karin, Z. Liu, E. Zandi, AP-1 function and regulation, Curr Opin Cell Biol 9(2) (1997) 240-6.
P.W. Vesely, P.B. Staber, G. Hoefler, L. Kenner, Translational regulation mechanisms of AP-1 proteins, Mutat Res 682(1) (2009) 7-12.
E. Shaulian, M. Karin, AP-1 in cell proliferation and survival, Oncogene 20(19) (2001) 2390-400.
M. Kappelmann, A. Bosserhoff, S. Kuphal, AP-1/c-Jun transcription factors: regulation and function in malignant melanoma, Eur J Cell Biol 93(1-2) (2014) 76-81.
J.J. Ventura, N.J. Kennedy, R.A. Flavell, R.J. Davis, JNK regulates autocrine expression of TGF-beta1, Mol Cell 15(2) (2004) 269-78.
A. Tokito, M. Jougasaki, Matrix Metalloproteinases in Non-Neoplastic Disorders, Int J Mol Sci 17(7) (2016).
N. Martin, O. Soriani, D. Bernard, Cardiac Glycosides as Senolytic Compounds, Trends Mol Med 26(3) (2020) 243-245.
S. He, N.E. Sharpless, Senescence in Health and Disease, Cell 169(6) (2017) 1000-1011.
M. Paez-Ribes, E. Gonzalez-Gualda, G.J. Doherty, D. Munoz-Espin, Targeting senescent cells in translational medicine, EMBO Mol Med 11(12) (2019) e10234.

Download PDF

Journal Publication

published 14 Apr, 2021

Read the published version in Journal of Orthopaedic Surgery and Research →

Editorial decision: Major revision
16 Nov, 2020
Review #3 received at journal
28 Sep, 2020
Review #1 received at journal
27 Sep, 2020
Reviewer #3 agreed at journal
14 Sep, 2020
Reviewer #2 agreed at journal
16 Aug, 2020
Reviewer #1 agreed at journal
14 Aug, 2020
Reviewers invited by journal
13 Aug, 2020
Editor assigned by journal
27 Jul, 2020
Editor invited by journal
26 Jul, 2020
Submission checks completed at journal
22 Jul, 2020
First submitted to journal
21 Jul, 2020

You are reading this older preprint version

Read the latest preprint version →

Identification of Key Regulators Responsible for Dysregulated Networks in Osteoarthritis by Large-Scale Expression Analysis

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Methods

Literature search and data collection

Data normalization and removal of batch effects

Filtering of differentially expressed genes

Integration of the PPI network and hub gene identification

Identification of putative target genes of JUN

Statistical analysis

Results

Outline of data analysis

Molecular characteristic of OA

Protein-Protein interaction network of DEGs

JUN functions as a key TF that associated with the development of OA

Identification and characterization of molecular subgroups

Drug repurposing for OA

Discussion

Conclusion

Abbreviations

Declarations

Acknowledgements

Author Contributions

Funding

Availability of data and materials

Ethics approval and consent to participate.

Consent for publication

Competing interests

References

Supplementary Files

Status:

Journal Publication

Version 1