BioMed Research International
Volume 2020 (2020), Article ID 2379280, 10 pages
https://doi.org/10.1155/2020/2379280

Research Article

Study of Osteoarthritis-Related Hub Genes Based on Bioinformatics Analysis

Zhengqing Zhu, Lei Zhong, Ronghang Li, Yuzhe Liu, Xiangrun Chen, Zhaoyan Li, and Lanfeng Huang

Departments of Orthopaedics, The Second Hospital of Jilin University, Changchun, Jilin, China

Correspondence should be addressed to Zhaoyan Li; zyli17@mails.jlu.edu.cn and Lanfeng Huang; hlf@jlu.edu.cn

Received 24 April 2020; Revised 2 July 2020; Accepted 17 July 2020; Published 6 August 2020

Academic Editor: Betti Giusti

Copyright © 2020 Zhengqing Zhu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Osteoarthritis (OA) is a common cause of morbidity and disability worldwide. However, the pathogenesis of OA is unclear. Therefore, this study was conducted to characterize the pathogenesis and implicated genes of OA. The gene expression profiles of GSE82107 and GSE55235 were downloaded from the Gene Expression Omnibus database. Altogether, 173 differentially expressed genes including 68 upregulated genes and 105 downregulated genes in patients with OA were selected based on the criteria of and an adjusted value < 0.05. Protein-protein interaction network analysis showed that FN1, COL1A1, IGF1, SPP1, TIMP1, BGN, COL5A1, MMP13, CLU, and SDC1 are the top ten genes most closely related to OA. Quantitative reverse transcription-polymerase chain reaction showed that the expression levels of COL1A1, COL5A1, TIMP1, MMP13, and SDC1 were significantly increased in OA. This study provides clues for the molecular mechanism and specific biomarkers of OA.

1. Introduction

Osteoarthritis (OA) is a common orthopedic disease that occurs in middle-aged and elderly people [1, 2]. OA is characterized by joint degeneration, with the main clinical manifestations being pain, swelling, and limited movement of the involved joints. Its etiology is related to joint injury, joint or limb dysplasia, age, obesity, and genetic factors. OA is among the greatest causes of morbidity and disability worldwide [3]. According to the estimation, it is estimated that approximately 31 million people in the United States are affected [4]. Local or systemic inflammation is considered to play a key role in the progress of OA, particularly in the early stage of disease [5]. Although abnormal gene expression in OA synovial cells has been widely reported in human and animal studies, altered gene expression and regulation in OA synovial cells remain unclear [6].

In recent years, a lot of research had been conducted to explore the molecular characteristics of OA. High-throughput microarray methods have attracted attention and are widely used in fields such as medical oncology ranging from diagnosis to treatment and prediction of the pathogenesis and prognosis of patients and for identifying drug targets [7–9]. Epigenetic regulation of gene expression in OA progression has recently been reported, and bioinformatics and microarray technologies have been extensively conducted to screen gene signaling pathways of OA [10, 11], providing strong theoretical support for the diagnosis as well as the treatment of OA [12, 13]. The roles of some differentially expressed genes (DEGs) in OA were studied by bioinformatics analysis and microscopy, enabling studies of the complex process of OA occurrence and development [14]. However, microarray methods are limited by the small number of samples collected, measurement errors, and inadequate information collection. The integrative analysis of multiple factors in OA remains challenging, and the exact pathogenesis of this condition remains unclear.

Despite these limitations, microarray systems are still efficient tools for detecting gene expression, identifying biomarkers, and evaluating epigenetic variation [15]. Network modeling of protein-protein interactions (PPIs) is a new technology for studying diseases and identifying disease-related gene targets [16]. To further explore DEGs and their molecular biological roles in OA, we have downloaded raw microarray data (GSE82107, GSE55235) from the Gene Expression Omnibus (GEO). DEGs between osteoarthritis synovium and normal controls were screened out by a bioinformatics method. A PPI network was built to enrich the protein domain of genes in the PPI network module and screen for hub genes related to OA. Significant DEGs between patients with OA and normal subjects may play an important part in OA development and progression. Additionally, the understanding of the underlying molecular mechanism of the pathogenesis of OA may be improved by using this method, revealing new therapeutic approaches for the epigenetic regulation of OA. In this research, we intended to explore the pathogenesis of OA through gene expression analysis to identify new potential biomarkers.

2. Materials and Methods

2.1. Microarray Data Source

The GEO dataset is a public functional genomics data repository that stores public gene expression datasets and platform records [17]. The gene expression profiles of GSE82107 and GSE55235 were downloaded from the GEO database. The GSE82107 dataset was obtained by using the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array, Santa Clara, CA, USA); the microarray data of GSE82107 contained 20 synovium samples from patients with OA and 17 knee synovium samples from normal subjects [18]. The GSE55235 dataset was based on the GPL96 platform (Affymetrix Human Genome U133A Array), and the microarray data of GSE55235 contained 10 OA patients synovium and 10 normal knee synovium [19].

2.2. Identification of Differentially Expressed Genes

Raw data were corrected using the RMA method before analyzing DEGs in OA, and we used a limma package in R software to screen for DEGs between the synovium of patients with OA and normal knee synovium [20]. DEGs were screened according to and an adjusted value < 0.05. Common DEGs in the GSE82107 and GSE55235 datasets were screened with FunRich software [21].

2.3. GO and Pathway Enrichment Analysis

To better understand the pathways and signal transduction processes involving DEGs in disease processes, the online DAVID bioinformatics database was used for GO and KEGG pathway analysis, and a gene and were considered to indicate significant results [22]. The GO database is the world’s largest source of information on gene function, with enrichment analysis including the molecular function (MF), cellular component (CC), and categories of biological processes (BP). [23].

2.4. PPI Network Construction

The STRING database is used to integrate all publicly available sources of PPI information and is also used to calculate predictions based on PPI analysis [24]. We utilized the STRING database to structure a PPI network of DEGs. An interaction was selected as a significant threshold. We led the raw data in Cytoscape software and used the cytoHubba plugin to build a subnetwork of PPIs and screened hub genes [25].

2.5. Case and Control Groups

Our research was approved by the ethics committee of the Second Hospital of Jilin University (Ethics number 2018-292). Ten patients with OA and ten patients with meniscus tears were included in this study, and all patients signed informed consent. We collected OA synovial tissue from patients undergoing total knee arthroplasty at the Jilin University Second Hospital, and normal synovial tissue was collected from patients undergoing arthroscopic surgery at the Jilin University Second Hospital.

2.6. Validation of Gene Expression

Real-time PCR was worked to validate the top ten hub genes, and total RNA was extracted from osteoarthritic and control synovial tissues by utilizing TRIzol reagent (Invitrogen, Carlsbad, CA, USA) followed by reverse transcription into cDNA. The mRNA expression levels were further assessed by quantitative reverse transcription- (qRT-) PCR using the QuantStudio™ 7 Flex real-time PCR system (Applied Biosystems, Foster City, CA, USA). These primers were designed with primer premier 6.0 software (Table 1). All samples were standardized to GAPDH expression; also, the experiment was repeated three times. The relative abundance of genes was calculated using the 2-ΔΔCt method, and data were analyzed with GraphPad software (GraphPad, Inc., La Jolla, CA, USA). The value < 0.05 was considered significant. The Pearson correlation coefficient was used to examine the relationship between key genes.

Table 1: The primers of the top 10 hub genes.

Gene
Forward primer
Reverse primer

FN1
AGGAAGCCGAGGTTTTAACTG
AGGACGCTCATAAGTGTCACC

COL1A1
GAGGGCCAAGACGAAGACATC
CAGATCACGTCATCGCACAAC

IGF1
GCTCTTCAGTTCGTGTGTGGA
GCCTCCTTAGATCACAGCTCC

SPP1
GAAGTTTCGCAGACCTGACAT
GTATGCACCATTCAACTCCTCG

TIMP1
CTTCTGCAATTCCGACCTCGT
ACGCTGGTATAAGGTGGTCTG

BGN
CAGTGGCTTTGAACCTGGAG
GGGAGGTCTTTGGGGATGC

COL5A1
GCCCGGATGTCGCTTACAG
AAATGCAGACGCAGGGTACAG

MMP13
ACTGAGAGGCTCCGAGAAATG
GAACCCCGCATCTTGGCTT

CLU
CCAATCAGGGAAGTAAGTACGTC
CTTGCGCTCTTCGTTTGTTTT

SDC1
CTGCCGCAAATTGTGGCTAC
TGAGCCGGAGAAGTTGTCAGA

GAPDH
CGGACCAATACGACCAAATCCG
AGCCACATCGCTCAGACACC

3. Results

3.1. Identification of DEGs in OA

We evaluated 10 patients with OA and 10 healthy control patients. After gene expression data processing and normalizing, we utilized the limma package in R software to screen for DEGs in each GEO dataset. The GSE5235 and GSE82107 datasets were standardized, and the results are shown in Figures 1 and 2. By using and adjusted value < 0.05 as thresholds, we determined 553 upregulated genes and 972 downregulated genes in GSE55235 and also 2811 upregulated genes and 639 downregulated genes in GSE82107. Based on overlapping results in the Venn diagrams (Figure 3), 173 genes were identified; 68 genes were upregulated, and the other 105 genes were downregulated.

Figure 1: Boxplots of GSE55235 data prior to and following normalization. The red box plots represent the normalized data, and the green box plots represent the data before normalization.

Figure 2: Boxplots of GSE82107 data prior to and following normalization. The red box plots represent the normalized data, and the green box plots represent the data before normalization.

(a)

(b)

Figure 3: The Venn analysis of upregulated and downregulated DEGs. (a) The Venn analysis of upregulated DEGs showed there are 68 DEGs coexpressed in two gene profiles. (b) The Venn analysis of downregulated DEGs presented there are 105 DEGs that were coexpressed in two gene profiles.

3.2. GO and KEGG Pathway Enrichment Analyses

We conducted enrichment analysis to assess the biological processes and pathways within OA (an interaction was selected as a significance). The GO enrichment analysis results were divided into three functional categories, including MF, CC, and BP. For BP, the DEGs were enriched in cell adhesion, immune response, and osteoblast differentiation (Table 2). In the CC group, the DEGs were majorly enriched in the extracellular space, extracellular matrix, and extracellular exosome. In the MF group, DEGs were majorly enriched in the collagen binding, integrin binding, and extracellular matrix structural constituent (Figure 4). We used the DAVID online database to perform DEG pathway enrichment analysis. KEGG analysis showed that DEGs play a role in regulating synovial inflammation in OA through various pathways including the AMPK signaling pathway, osteoclast differentiation, insulin signaling pathway, autophagy, ECM-receptor interaction, and HIF-1 signaling pathway (Figure 5, Table 3).

Table 2: The significantly enriched analysis of DEGs in OA.

Expression
Category
Term
Description
Gene count
value

DEGs
BP
BP
BP
BP
BP
CC
CC
CC
CC
CC
MF
MF
MF
MF
MF
GO:0007155
GO:0030198
GO:0050870
GO:0001525
GO:0001649
GO:0005615
GO:0031012
GO:0070062
GO:0005576
GO:0005578
GO:0005518
GO:0005178
GO:0005201
GO:0016504
GO:0000977
Cell adhesion
Extracellular matrix organization
Positive regulation of T cell activation
Angiogenesis
Osteoblast differentiation
Extracellular space
Extracellular matrix
Extracellular exosome
Extracellular region
Proteinaceous extracellular matrix
Collagen binding
Integrin binding
Extracellular matrix structural constituent
Peptidase activator activity
RNA polymerase II regulatory region sequence-specific DNA binding
15
10
4
9
6
35
16
47
31
11
6
6
5
3
7
5.17-05
5.76-05
4.87-04
7.83-04
0.00220487
3.51-09
2.52-08
2.26-06
1.95-05
8.13-05
1.79-04
0.002314609
0.002834025
0.003280773
0.010113065

Figure 4: According to gene ontology analysis, DEGs can be divided into three categories: molecular function, biological process, and cellular component.

Figure 5: Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the pathways. The size of the black spots stands for the gene number; the gradual color stands for the value.

Table 3: Signaling pathway enrichment analysis of DEG function in OA.

Expression
Term
Description
Gene count
value

DEGs
hsa04512
ECM-receptor interaction
6
0.001503576

hsa04066
HIF-1 signaling pathway
5
0.013969316

hsa04510
Focal adhesion
7
0.01467336

hsa04152
AMPK signaling pathway
5
0.031444563

hsa04380
Osteoclast differentiation
5
0.038325066

hsa04910
Insulin signaling pathway
5
0.045003644

hsa04145
Phagosome
5
0.057882463

hsa05202
Transcriptional misregulation in cancer
5
0.079165687

hsa01212
Fatty acid metabolism
3
0.079189222

3.3. Identification of Hub Genes

To further examine the hub genes involved in the development of OA, we used the STRING database to evaluate the interaction between genes. We built a PPI network by using Cytoscape software and some data from the STRING database. To further analyze protein interactions, we evaluated the betweenness centrality and degree using the plugin cytoHubba in Cytoscape software (Figure 6). Hub genes with betweenness centrality and degree indicate that these genes play a key role in this network. The 10 hub genes showing significant interactions were FN1, COL1A1, IGF1, SPP1, TIMP1, BGN, COL5A1, MMP13, CLU, and SDC1 (Figure 7).

Figure 6: PPI network constructed with the downregulated and upregulated DEGs. Blue nodes represent upregulated genes, and red nodes represent downregulated genes. The size of the triangles and circles represents the degree value.

Figure 7: The 10 hub genes showed significant interactions, and these genes were validated by qRT-PCR. (The color of the boxes is the weight of hub genes. From red to yellow, red indicates a higher degree of association.)

3.4. Validation of Hub Genes

To check the data analysis results, qRT-PCR was used to detect the expression level of the first 10 hub genes in the OA synovium of the knee joint and normal control groups. Statistical analysis proved that ADIPOQ, IL6, and CXCR1 were obviously raised in the synovium of OA samples () (Figure 8). All validations are consistent with the results of this research.

Figure 8: Validation of the top10 hub genes by qRT-PCR between the OA group () and the control group (). The dotted line at ordinate 1 represented the relative mRNA expression of the control group. All samples were normalized to the expression of GAPDH, and the relative expression levels of each gene were analyzed using the 2−ΔΔCt method. , , and .

4. Discussion

OA is a complex disease caused by genetic and environmental factors. This condition affects more than 40% of 70 year olds and is the main cause of loss of body movement and pain [26]. OA is characterized by several factors, including the loss of articular cartilage, hyperosteogeny, and synovitis. The combined effects of mechanical and biochemical factors are implicated in OA development. The rate of OA increases with age and mechanical wear on the joints. However, the exact pathogenesis of OA is unclear, and there are no effective treatment methods; joint replacement surgery is the last resort for treating OA [27].

Gene chip and high-throughput sequencing technologies can be used for detecting the gene expression, microRNA, long noncoding RNA, and DNA methylation to explore genetic alterations in disease [28]. Microarray techniques have also been widely used to predict potential target genes for OA [7, 29]. To decrease the number of false-positive results, we performed functional enrichment and network analysis of the DEGs. A total of 173 genes were identified, with 68 genes upregulated and 105 genes downregulated.

Using the KEGG signal pathway analysis method to analyze the enrichment of DEGs. Based on the string online database, a PPI network is built by using the software of Cytoscape, and 10 most closely related genes were screened out by the degree analysis method of the cytoHubba plugin. The following genes were highly expressed: FN1, COL1A1, IGF1, SPP1, TIMP1, BGN, COL5A1, MMP13, CLU, and SDC1. To validate these findings, the qRT-PCR results confirmed that the expression levels of COL1A1, COL5A1, TIMP1, MMP13, and SDC1 were obviously increased in OA samples ().

The main organic component of bone is collagen type I, which consists of one pro-alpha2 (I) chain encoded by the collagen type I alpha 2 chain gene (COL1A2) and two pro-alpha1 (I) chains encoded by two the collagen type I alpha 1 chain genes (COL1A1) [30]. COL1A1 is one of several collagens showing rich differences in RNA and protein levels. Collagens are the major structural components of the cartilage. Collagen imbalance plays an important pathogenic role in OA [31, 32]. Type I collagen is the most abundant collagen in scar tissue and is the final product of tissue healing. A recent study showed that COL1A1 was upregulated in the synovium of patients with advanced OA as well as in the synovium of OA-induced mice and human fibroblasts stimulated by transforming growth factor-β (TGF-β) [33]. According to previous studies, collagen type I is rare in normal chondrocytes but is highly expressed in OA cartilage, especially in the late stage of OA [34]. The collagen type I alpha 5 chain gene (COL5A1) encodes the α1 chain of type V collagen and is small fibrous collagen found in ligaments, tendons, and other tissues [35]. Type V collagen, encoded by COL5A1, is a type I collagen present in tissues containing type I collagen and may play an important role in regulating the assembly of heteromorphic fibers composed of type I as well as type V collagen [36, 37]. Because these genes were found to be upregulated in OA patients, collagen may play a key role in the pathogenesis of OA.

It has been established that the aging of chondrocytes in elderly patients promotes the development of OA [38]. In this case, aging and pathological cells still survive but show changes in the matrix metalloproteinase (MMP) expression profile. MMPs are known to be involved in cartilage destruction and inflammation in the joint [39]. MMPs belong to a protease family called metzincin superfamily, which contains 23 core members and more than 30 extended members [40]. These enzymes participated in all kinds of physiological and pathological processes in life activities, such as morphogenesis, tissue repair, wound healing, and postinjury remodeling, as well as the progression of OA, chronic tissue ulcers, cancer, and arthritis [41]. MMP13 is a member of the MMP family. Due to its high expression in human OA cartilage and the ability to degrade type II collagen fibers, MMP13 is supposed to the main collagenase involved in the growth of OA [42–44]. In the past two decades, studies have shown that MMP13 is significant in the degradation of cartilage aggregation proteins and is closely related to the senescence of chondrocytes in the pathogenesis of OA [45, 46]. Therefore, MMP13 is considered a promising target for OA therapy [47]. The tissue inhibitor of metalloproteinase (TIMP) proteins forms a family of natural inhibitors of MMPs. This family consists of four members: TIMP 1, TIMP 2, TIMP3, and TIMP 4 [48, 49]. Because of the strong binding affinity between TIMP1 and MMP13, TIMP1 is a functional regulator of MMP13-related chondrocyte aging and cartilage matrix changes. Therefore, the improvement of TIMP1 activity is considered to be an effective way to treat OA [50].

Syndecan 1 (SDC1) protein is a transmembrane (type I) heparin sulfate proteoglycan, which is a member of the syndecan proteoglycan family [51]. SDC1 protein has several molecular roles in proliferation, migration, and cell-matrix interactions through its extracellular matrix protein receptor. In the early stage of cartilage degeneration in OA, the expression of SDC1 in the articular cartilage was upregulated, indicating that SDC1 participates in the repair of cartilage fibrils [51, 52].

In conclusion, the PPI network and data analyses are effective ways for screening OA target proteins and corresponding gene targets, which can further reveal the pathological and biological mechanisms of OA. In the present study, based on qRT-PCR results, COL1A1, COL5A1, TIMP1, MMP13, and SDC1 may be gene diagnostic markers and drug targets of OA. Even though we validated several hub genes and pathways, the small sample size is one limitation of this study. Further large-scale experimental studies are needed to confirm these preliminary results.

Data Availability

All data generated or analyzed during this study are included in this article.

Conflicts of Interest

The authors declare that there is no conflict of interests regarding the publication of this article.

Authors’ Contributions

All authors were involved in the study and approved the final version of the manuscript.

Acknowledgments

This work was supported by the Health Special Project of Jilin Province (Grant Number 201817294302), the Department of Science and Technology of Jilin Province (Grant Number 20200201327JC), the Education Department of Jilin Province (Grant Number 3D518Q243429), and the Training Program of Outstanding Doctoral Student by Norman Bethune Health Science Center of Jilin University (Grant Number 470110000646). We thank Dr. Zhaoyan Li for the assistance with the samples.

References

R. Higgs, “Concentrated efforts to detect early OA,” Nature Reviews Rheumatology, vol. 6, no. 11, p. 616, 2010.

D. G. Hoy, E. Smith, M. Cross et al., “The global burden of musculoskeletal conditions for 2010: an overview of methods,” Annals of the Rheumatic Diseases, vol. 73, no. 6, pp. 982–989, 2014.

L. Parkinson, D. L. Waters, and L. Franck, “Systematic review of the impact of osteoarthritis on health outcomes for comorbid disease in older people,” Osteoarthritis and Cartilage, vol. 25, no. 11, pp. 1751–1770, 2017.

J. Park, A. Mendy, and E. R. Vieira, “Various types of arthritis in the United States: prevalence and age-related trends from 1999 to 2014,” American Journal of Public Health, vol. 108, no. 2, pp. 256–258, 2018.

W. H. Robinson, C. M. Lepus, Q. Wang et al., “Low-grade inflammation as a key mediator of the pathogenesis of osteoarthritis,” Rheumatology, vol. 12, no. 10, pp. 580–592, 2016.

M. Zhang, K. Lygrissea, and J. Wanga, “Role of microRNA in osteoarthritis,” Journal of Arthritis, vol. 6, no. 2, 2017.

Z. Li, Q. Wang, G. Chen et al., “Integration of gene expression profile data to screen and verify hub genes involved in osteoarthritis,” BioMed Research International, vol. 2018, Article ID 9482726, 10 pages, 2018.

Z. Li, L. Zhong, Z. Du et al., “Network analyses of differentially expressed genes in osteoarthritis to identify hub genes,” BioMed Research International, vol. 2019, Article ID 8340573, 9 pages, 2019.

M. Nannini, M. A. Pantaleo, A. Maleddu, A. Astolfi, S. Formica, and G. Biasco, “Gene expression profiling in colorectal cancer using microarray technologies: results and perspectives,” Cancer Treatment Reviews, vol. 35, no. 3, pp. 201–209, 2009.

The Gene Ontology Consortium, “Expansion of the Gene Ontology knowledgebase and resources,” Nucleic Acids Research, vol. 45, no. D1, pp. D331–D338, 2017.

F. Q. Liu, “Analysis of differentially expressed genes in rheumatoid arthritis and osteoarthritis by integrated microarray analysis,” Journal of Cellular Biochemistry, vol. 120, no. 8, pp. 12653–12664, 2019.

M. Kanehisa, M. Furumichi, M. Tanabe, Y. Sato, and K. Morishima, “KEGG: new perspectives on genomes, pathways, diseases and drugs,” Nucleic Acids Research, vol. 45, no. D1, pp. D353–D361, 2017.

W. C. Li, D. L. Bai, Y. Xu et al., “Identification of differentially expressed genes in synovial tissue of rheumatoid arthritis and osteoarthritis in patients,” Journal of Cellular Biochemistry, vol. 120, no. 3, pp. 4533–4544, 2019.

M. C. Hochberg, L. Yerges-Armstrong, M. Yau, and B. D. Mitchell, “Genetic epidemiology of osteoarthritis: recent developments and future directions,” Current Opinion in Rheumatology, vol. 25, no. 2, pp. 192–197, 2013.

J. Liu, H. Li, L. Sun, Z. Wang, C. Xing, and Y. Yuan, “Aberrantly methylated-differentially expressed genes and pathways in colorectal cancer,” Cancer Cell International, vol. 17, no. 1, p. 75, 2017.

Y. Liu, S. Patel, R. Nibbe et al., “Systems biology analyses of gene expression and genome wide association study data in obstructive sleep apnea,” Pacific Symposium on Biocomputing, pp. 14–25, 2011.

T. Barrett, S. E. Wilhite, P. Ledoux et al., “NCBI GEO: archive for functional genomics data sets--update,” Nucleic Acids Research, vol. 41, Database issue, pp. D991–D995, 2013.

M. G. A. Broeren, M. de Vries, M. B. Bennink et al., “Functional tissue analysis reveals successful cryopreservation of human osteoarthritic synovium,” PLoS One, vol. 11, no. 11, article e0167076, 2016.

D. Woetzel, R. Huber, P. Kupfer et al., “Identification of rheumatoid arthritis and osteoarthritis patients by transcriptome-based rule set generation,” Arthritis Research & Therapy, vol. 16, no. 2, p. R84, 2014.

M. E. Ritchie, B. Phipson, D. Wu et al., “Limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic Acids Research, vol. 43, no. 7, article e47, 2015.

M. Pathan, S. Keerthikumar, C. S. Ang et al., “FunRich: an open access standalone functional enrichment and interaction network analysis tool,” Proteomics, vol. 15, no. 15, pp. 2597–2601, 2015.

D. W. Huang, B. T. Sherman, Q. Tan et al., “The DAVID Gene Functional Classification Tool: a novel biological module-centric algorithm to functionally analyze large gene lists,” Genome Biology, vol. 8, no. 9, article R183, 2007.

M. Ashburner, C. A. Ball, J. A. Blake et al., “Gene ontology: tool for the unification of biology. The Gene Ontology Consortium,” Nature Genetics, vol. 25, no. 1, pp. 25–29, 2000.

D. Szklarczyk, A. L. Gable, D. Lyon et al., “STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets,” Nucleic Acids Research, vol. 47, no. D1, pp. D607–D613, 2019.

C. H. Chin, S. H. Chen, H. H. Wu, C. W. Ho, M. T. Ko, and C. Y. Lin, “cytoHubba: identifying hub objects and sub-networks from complex interactome,” BMC Systems Biology, vol. 8, article S11, Supplement 4, 2014.

T. Vos, A. D. Flaxman, M. Naghavi et al., “Years lived with disability (YLDs) for 1160 sequelae of 289 diseases and injuries 1990-2010: a systematic analysis for the Global Burden of Disease Study 2010,” The Lancet, vol. 380, no. 9859, pp. 2163–2196, 2012.

P. A. Dieppe and L. S. Lohmander, “Pathogenesis and management of pain in osteoarthritis,” Lancet, vol. 365, no. 9463, pp. 965–973, 2005.

Z. Dong, J. Wang, T. Zhan, and S. Xu, “Identification of prognostic risk factors for esophageal adenocarcinoma using bioinformatics analysis,” Oncotargets and Therapy, vol. Volume 11, pp. 4327–4337, 2018.

Y. Fang, P. Wang, L. Xia et al., “Aberrantly hydroxymethylated differentially expressed genes and the associated protein pathways in osteoarthritis,” PeerJ, vol. 7, article e6425, 2019.

H. Zhang, C. Chen, Y. Cui et al., “lnc-SAMD14-4 can regulate expression of the COL1A1 and COL1A2 in human chondrocytes,” PeerJ, vol. 7, article e7491, 2019.

E. V. Tchetina, “Developmental mechanisms in articular cartilage degradation in osteoarthritis,” Arthritis, vol. 2011, Article ID 683970, 16 pages, 2011.

B. Xia, di Chen, J. Zhang, S. Hu, H. Jin, and P. Tong, “Osteoarthritis pathogenesis: a review of molecular mechanisms,” Calcified Tissue International, vol. 95, no. 6, pp. 495–505, 2014.

D. F. G. Remst, A. B. Blom, E. L. Vitters et al., “Gene expression analysis of murine and human osteoarthritis synovium reveals elevation of transforming growth factor β-responsive genes in osteoarthritis-related fibrosis,” Arthritis & Rheumatology, vol. 66, no. 3, pp. 647–656, 2014.

N. Miosge, M. Hartmann, C. Maelicke, and R. Herken, “Expression of collagen type I and type II in consecutive stages of human osteoarthritis,” Histochemistry and Cell Biology, vol. 122, no. 3, pp. 229–236, 2004.

D. E. Birk, “Type V collagen: heterotypic type I/V collagen interactions in the regulation of fibril assembly,” Micron : the international research and review journal for microscopy, vol. 32, no. 3, pp. 223–237, 2001.

M. Colombi, C. Dordoni, M. Venturini et al., “Spectrum of mucocutaneous, ocular and facial features and delineation of novel presentations in 62 classical Ehlers-Danlos syndrome patients,” Clinical Genetics, vol. 92, no. 6, pp. 624–631, 2017.

K. Willard, S. Mannion, C. J. Saunders, M. Collins, and A. V. September, “The interaction of polymorphisms in extracellular matrix genes and underlying miRNA motifs that modulate susceptibility to anterior cruciate ligament rupture,” Journal of Science and Medicine in Sport, vol. 21, no. 1, pp. 22–28, 2018.

J. S. Price, J. G. Waters, C. Darrah et al., “The role of chondrocyte senescence in osteoarthritis,” Aging Cell, vol. 1, no. 1, pp. 57–65, 2002.

B. J. Rose and D. L. Kooyman, “A tale of two joints: the role of matrix metalloproteases in cartilage biology,” Disease Markers, vol. 2016, Article ID 4895050, 7 pages, 2016.

H. Nagase, R. Visse, and G. Murphy, “Structure and function of matrix metalloproteinases and TIMPs,” Cardiovascular Research, vol. 69, no. 3, pp. 562–573, 2006.

G. Murphy and H. Nagase, “Progress in matrix metalloproteinase research,” Molecular Aspects of Medicine, vol. 29, no. 5, pp. 290–308, 2008.

P. G. Mitchell, H. A. Magna, L. M. Reeves et al., “Cloning, expression, and type II collagenolytic activity of matrix metalloproteinase-13 from human osteoarthritic cartilage,” The Journal of Clinical Investigation, vol. 97, no. 3, pp. 761–768, 1996.

A. Tsezou, “Osteoarthritis year in review 2014: genetics and genomics,” Osteoarthritis and Cartilage, vol. 22, no. 12, pp. 2017–2024, 2014.

K. Yamamoto, H. Okano, W. Miyagawa et al., “MMP-13 is constitutively produced in human chondrocytes and co-endocytosed with ADAMTS-5 and TIMP-3 by the endocytic receptor LRP1,” Matrix biology : journal of the International Society for Matrix Biology, vol. 56, pp. 57–73, 2016.

C. B. Forsyth, A. Cole, G. Murphy, J. L. Bienias, H. J. Im, and R. F. Loeser, “Increased matrix metalloproteinase-13 production with aging by human articular chondrocytes in response to catabolic stimuli,” The Journals of Gerontology Series A: Biological Sciences and Medical Sciences, vol. 60, no. 9, pp. 1118–1124, 2005.

M. Inada, Y. Wang, M. H. Byrne et al., “Critical roles for collagenase-3 (Mmp13) in development of growth plate cartilage and in endochondral ossification,” Proceedings of the National Academy of Sciences of the United States of America, vol. 101, no. 49, pp. 17192–17197, 2004.

N. G. Li, Y. P. Tang, and J. A. Duan, “New selective inhibitors of MMP-13 for inflammatory diseases: a patent evaluation (W02012151158),” Expert Opinion on Therapeutic Patents, vol. 23, no. 5, pp. 669–675, 2013.

K. Brew and H. Nagase, “The tissue inhibitors of metalloproteinases (TIMPs): an ancient family with structural and functional diversity,” Biochimica et Biophysica Acta, vol. 1803, no. 1, pp. 55–71, 2010.

W. Zhang, C. Zhang, C. Luo, Y. Zhan, and B. Zhong, “Design, cyclization, and optimization of MMP13-TIMP1 interaction-derived self-inhibitory peptides against chondrocyte senescence in osteoarthritis,” International Journal of Biological Macromolecules, vol. 121, pp. 921–929, 2019.

V. Knäuper, S. Cowell, B. Smith et al., “The role of the C-terminal domain of human collagenase-3 (MMP-13) in the activation of procollagenase-3, substrate specificity, and tissue inhibitor of metalloproteinase interaction,” The Journal of Biological Chemistry, vol. 272, no. 12, pp. 7608–7616, 1997.

P. E. Barre, F. Redini, K. Boumediene, C. Vielpeau, and J. P. Pujol, “Semiquantitative reverse transcription-polymerase chain reaction analysis of syndecan-1 and -4 messages in cartilage and cultured chondrocytes from osteoarthritic joints,” Osteoarthritis and Cartilage, vol. 8, no. 1, pp. 34–43, 2000.

H. Salminen-Mankonen, A. M. Säämänen, M. Jalkanen, E. Vuorio, and L. Pirilä, “Syndecan-1 expression is upregulated in degenerating articular cartilage in a transgenic mouse model for osteoarthritis,” Scandinavian Journal of Rheumatology, vol. 34, no. 6, pp. 469–474, 2009.