BioMed Research International
Volume 2019 (2019), Article ID 8340573, 9 pages
https://doi.org/10.1155/2019/8340573

Research Article

Network Analyses of Differentially Expressed Genes in Osteoarthritis to Identify Hub Genes

Zhaoyan Li,1,2,3 Lei Zhong,1,3 Zhenwu Du,1,2,3 Gaoyang Chen,1,2,3 Jing Shang,1,3 Qiwei Yang,2,3 Qingyu Wang,1,2,3 Yang Song,1,3 and Guizhen Zhang1,2,3

¹Department of Orthopedics of the Second Hospital of Jilin University, Ziqiang Street 218, Changchun, Jilin 130041, China
²Research Centre of the Second Hospital of Jilin University, Ziqiang Street 218, Changchun, Jilin 130041, China
³The Engineering Research Centre of Molecular Diagnosis and Cell Treatment for Metabolic Bone Diseases of Jilin Province, Ziqiang Street 218, Changchun, Jilin 130041, China

Correspondence should be addressed to Yang Song; songyangjlu@jlu.edu.cn and Guizhen Zhang; zhangguizhenjlu@163.com

Received 11 March 2019; Revised 19 March 2019; Accepted 8 April 2019; Published 16 April 2019

Academic Editor: Magali Cucchiarini

Copyright © 2019 Zhaoyan Li 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.

Background. Osteoarthritis (OA) is the most common degenerative disease in orthopedics. However, the cause and underlying molecular mechanisms are not clear. This study aims to identify the hub genes and pathways involved in the occurrence of osteoarthritis. Methods. The raw data of GSE89408 were downloaded from the Gene Expression Omnibus (GEO) database, and the differentially expressed genes (DEGs) were identified by R software. The DAVID database was used for pathway and gene ontology analysis, and p<0.05 and gene count >2 were set as the cut-off point. Moreover, protein-protein interaction (PPI) network construction was applied for exploring the hub genes in osteoarthritis. The expression levels of the top ten hub genes in knee osteoarthritis synovial membranes and controls were detected by quantitative real-time PCR system. Results. A total of 229 DEGs were identified in osteoarthritis synovial membranes compared with normal synovial membranes, including 145 upregulated and 84 downregulated differentially expressed genes. The KEGG pathway analysis results showed that up-DEGs were enriched in proteoglycans in cytokine-cytokine receptor interaction, chemokine signaling pathway, rheumatoid arthritis, and TNF signaling pathway, whereas down-DEGs were enriched in the PPAR signaling pathway and AMPK signaling pathway. The qRT-PCR results showed that the expression levels of ADIPOQ, IL6, and CXCR1 in the synovium of osteoarthritis were significantly increased (p <0.05).

1. Introduction

Osteoarthritis is the most common degenerative disease worldwide; the global incidence of knee osteoarthritis is estimated at 3.8 percent in 2010 [1]. Osteoarthritis is the leading cause of disability in the US, with 43.5% of the arthritis patients having activity limitations [2]. Pain and stiffness, especially post-exercise pain and stiffness, are the main symptoms and have a considerable impact on the ability of daily life activities[1]. These symptoms severely reduce the quality of life and create a heavy socioeconomic burden.

Recent research has increasingly recognized the important role of the synovial membrane in the progression of osteoarthritis [3]. Synovitis in osteoarthritis is associated with activity disorder, joint severe pain, and cartilage loss in certain patient populations [4]. Osteoarthritis is a polygenic disorder. In recent years, genetic analysis data have revealed that OA-related genes are often associated with the development of synovial joints [5]. But until now, the hub gene and pathway of osteoarthritis are still not clear. Therefore, it is critical to elucidate the pathogenesis and progression of osteoarthritis further.

In this work, we have downloaded the microarray dataset GSE89408 from the GEO database, which contained 22 knee osteoarthritis synovial membranes and 28 matched normal synovial membranes. DEGs in osteoarthritis and normal synovial membranes were identified using the limma package in RStudio software, and DAVID database was used for pathway and gene ontology analysis of DEGs. Then, the PPI network was constructed to show the connection of DEGs and screening of hub genes. Our study aims to explore the hub genes and pathways of osteoarthritis, which may help to better understand the molecular mechanism of osteoarthritis.

2. Materials and Methods

2.1. Microarray Data Information

The GEO database is a public genomics data repository which provides high-throughput gene expression data [6]. The raw data of GSE89408 (GPL11154 platform, Illumina HiSeq 2000) was downloaded from the GEO database which includes 22 synovial membranes of knee osteoarthritis and 28 normal synovial membranes of controls [7].

2.2. Data Preprocessing and DEG Screening

The limma package in R was used to identify the DEGs between osteoarthritis synovial membranes and normal synovial membranes [8]. Background adjustment, normalization, and summarization were all included in the process of preprocessing. P values were adjusted using the Hochberg and Benjamini test, and p<0.05 and gene count >2 were set as the cut-off point.

2.3. Gene Ontology and Pathway Analysis

The DAVID Gene Functional Classification Tool can classify large genes into biological modules [9]. Gene Ontology provides the logical structure of the biological functions and their relationships to one another genes. GO analysis includes categories of molecular function (MF), cellular component (CC), and biological processes (BP). Pathway analysis is the process of classifying large genes by the KEGG database. In our study, candidate DEGs functions and pathway enrichment were analyzed using the DAVID database. P<0.05 and gene count >2 were set as the cut-off point.

2.4. PPI Network Integration

The STRING database is an online resource whose main function is to construct functional protein association networks [10]. The interaction scores >0.40 (medium confidence) were defined as significant. The PPI was constructed by Cytoscape software, and the interaction relationship of DEGs was analyzed by a plug-in of Cytoscape software. CytoHubba, as a plug-in Cytoscape software, provides an effective method for identifying hub genes in PPI networks through degree method [11].

2.5. qRT-PCR Validation and Statistical Analysis

Hub genes were verified by the qRT-PCR system. Total RNA was extracted from the control and knee osteoarthritis synovial membranes using TRIzol reagent (TaKaRa, Japan) and then reverse-transcribed to cDNA. Primers were designed by Primer 5.0 software (PREMIER Biosoft, Palo Alto, CA, USA), and a QuantStudio™ 7 Flex real-time PCR system (Applied Biosystems, Carlsbad, CA, USA) was used. Primers for mRNA are shown in Table 1. All samples were normalized to GAPDH. And the relative expression levels of each gene were calculated using 2−ΔΔCt methods. SPSS software (version 22.0 SPSS Inc.) was utilized to analyze statistical data, and P values < 0.05 were considered as statistically significant.

Table 1: The primers of top 10 hub genes.

Gene
Forward primer
Reverse primer

GAPDH
CGGACCAATACGACCAAATCCG
AGCCACATCGCTCAGACACC

IL8
ACTGAGAGTGATTGAGAGTGGAC
AACCCTCTGCACCCAGTTTTC

Il6
TCAATATTAGAGTCTCAACCCCCA
GAAGGCGCTTGTGGAGAAGG

ADIPOQ
GCTGGGAGCTGTTCTACTG
TACTCCGGTTTCACCGATGTC

CXCR1
CTGACCCAGAAGCGTCACTTG
CCAGGACCTCATAGCAAACTG

CXCL1
AACCGAAGTCATAGCCACAC
GTTGGATTTGTCACTGTTCAGC

LIPE
TCAGTGTCTAGGTCAGACTGG
AGGCTTCTGTTGGGTATTGGA

SLC2A4
TGGGCGGCATGATTTCCTC
GCCAGGACATTGTTGACCAG

FASN
AAGGACCTGTCTAGGTTTGATGC
TGGCTTCATAGGTGACTTCCA

FABP4
ACTGGGCCAGGAATTTGACG
CTCGTGGAAGTGACGCCTT

FPR1
TGGGAGGACATTGGCCTTTC
GGATGCAGGACGCAAACAC

2.6. Patients and Controls

Our research was approved by the Ethics Committee of the Second Hospital of Jilin University, Jilin, China. 10 knee osteoarthritis patients and 10 meniscal tear patients without obvious synovitis were enrolled, and all gave informed consent. Synovial samples of knee osteoarthritis were obtained from patients undergoing total knee arthroplasty in the Second Hospital of Jilin University, Jilin, China. And normal synovial membranes were obtained from arthroscopic surgery cases upon meniscus injury in the Second Hospital of Jilin University, Jilin, China.

3. Results

3.1. Differentially Expressed Genes

The GSE89408 dataset was standardized, and the results are shown in Figure 1. A total of 22 knee osteoarthritis synovial membranes and 28 matched normal samples were analyzed, using and p < 0.05 as the cut-off criteria, and we ultimately obtained 229 DEGs in osteoarthritis synovial membranes compared with normal synovial membranes, including 145 upregulated and 84 downregulated DEGs.

Figure 1: Gene expression before and after normalization: The green box plots represent the data before normalization, and the red box plots represent the normalized data.

3.2. Enrichment Analysis of DEGs

To acquire the functions of DEGs, GO function enrichment was analyzed by DAVID database. The top five BP terms, CC terms, and MF terms are shown in Table 2. In the BP group, the up-DEGs are mainly enriched in the inflammatory response, cell-cell signaling, chemokine-mediated signaling pathway, and immune response, and the down-DEGs are mainly enriched in lipid metabolic process and brown fat cell differentiation. In the CC group, the up-DEGs are mainly enriched in the extracellular region and extracellular space, and the down-DEGs are mainly enriched in lipid particle and mitochondrion. Moreover, in the MF group, the up-DEGs are mainly enriched in CXCR chemokine receptor binding, RAGE receptor binding, cytokine activity, and chemokine activity, and the down-DEGs are mainly enriched in drug binding and protein homodimerization activity. The results of GO enrichment analysis are shown in Figures 2 and 3. The pathway analysis results showed that up-DEGs were enriched in proteoglycans in cytokine-cytokine receptor interaction, chemokine signaling pathway, rheumatoid arthritis, and TNF signaling pathway, whereas down-DEGs were enriched in the PPAR signaling pathway and AMPK signaling pathway. The results are shown in Figure 4 and Table 3.

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

Expression
Category
Term
Description
Gene Count
P Value

UP-DEGs
BP
GO:0006954
inflammatory response
19
1.21E-13

BP
GO:0007267
cell-cell signaling
11
3.03E-07

BP
GO:0070098
chemokine-mediated signaling pathway
6
2.07E-05

BP
GO:0006935
chemotaxis
7
2.34E-05

BP
GO:0006955
immune response
11
2.71E-05

CC
GO:0005615
extracellular space
25
1.66E-09

CC
GO:0005576
extracellular region
27
2.28E-09

CC
GO:0005578
proteinaceous extracellular matrix
6
0.006499621

CC
GO:0016323
basolateral plasma membrane
4
0.044790466

CC
GO:0031988
membrane-bounded vesicle
2
0.05975494

MF
GO:0008009
chemokine activity
5
6.54E-05

MF
GO:0045236
CXCR chemokine receptor binding
3
6.87E-04

MF
GO:0050786
RAGE receptor binding
3
0.001043791

MF
GO:0005540
hyaluronic acid binding
3
0.004638666

MF
GO:0005125
cytokine activity
5
0.007811921

DOWN-DEGs
BP
GO:0006629
lipid metabolic process
5
0.001123386

BP
GO:0050873
brown fat cell differentiation
3
0.003912643

BP
GO:0010890
positive regulation of sequestering of triglyceride
2
0.020252031

BP
GO:0010642
negative regulation of platelet-derived
2
0.020252031

growth factor receptor signaling pathway

BP
GO:0042593
glucose homeostasis
3
0.035075126

CC
GO:0005811
lipid particle
7
6.93E-08

CC
GO:0005739
mitochondrion
11
0.009807547

CC
GO:0070062
extracellular exosome
16
0.031369825

CC
GO:0005576
extracellular region
11
0.032718894

CC
GO:0005615
extracellular space
9
0.067564494

MF
GO:0042803
protein homodimerization activity
7
0.016759119

MF
GO:0008144
drug binding
3
0.01974095

Table 3: Signaling pathway enrichment analysis of DEGs function in osteoarthritis.

Expression
Term
Description
Gene Count
P Value

DOWN-DEGs
hsa04923
Regulation of lipolysis in adipocytes
4
7.05E-04

hsa03320
PPAR signaling pathway
4
0.001191375

hsa04152
AMPK signaling pathway
4
0.006702552

hsa04910
Insulin signaling pathway
3
0.07107271

UP-DEGs
hsa04060
Cytokine-cytokine receptor interaction
8
1.22E-04

hsa04062
Chemokine signaling pathway
6
0.001756326

hsa05202
Transcriptional misregulation in cancer
5
0.007918451

hsa05323
Rheumatoid arthritis
4
0.008366703

hsa04668
TNF signaling pathway
4
0.014229192

hsa05322
Systemic lupus erythematosus
4
0.02575911

hsa05150
Staphylococcus aureus infection
3
0.027341561

Figure 2: The top ten GO terms in each group: The top ten GO terms in each group were ranked by p value. The gradual color represents the z-score.

Figure 3: GO enrichment analysis of DEGs: The gradual color represents the . The genes were ordered according to their values setting gene.

Figure 4: KEGG enrichment analysis of the pathways: The gradual color represents the P value, the size of the black spots represents the gene number.

3.3. PPI Network Analysis and Hub Genes Screening

Using the STRING database and Cytoscape software, we get a total of 67 nodes including 42 upregulated and 25 downregulated DEGs (Figure 5). Through plug-in CytoHubba in Cytoscape software, we evaluated the degree and betweenness centrality in the PPI network and screening the hub genes. The 10 hub genes showing significant interaction were IL8, ADIPOQ, CXCR1, CXCL1, LIPE, FPR1, FABP4, SLC2A4, FASN, and IL6 (Figure 6).

Figure 5: PPI network of DEGs: Red circles represent upregulated DEGs, and blue squares represent downregulated DEGs. The size of the nodes represents the degree value.

Figure 6: PPI network of hub genes: The gradual color represents the degree value. The circles represent upregulated genes, and squares represent downregulated genes.

3.4. Validation of Hub Genes

To verify the results of microarray, the expression levels of the top ten hub genes in knee osteoarthritis synovial membranes and controls were detected by the qRT-PCR system. The statistical results showed that the expression levels of ADIPOQ, IL6, and CXCR1 in the synovium of osteoarthritis were significantly increased (p <0.05) (Figure 7). All validations are consistent with the analytical results in this study.

Figure 7: Validation of the top 10 hub genes by qRT-PCR between the OA group and the control group: All samples were normalized to GAPDH. And the relative expression levels of each gene were calculated using 2−ΔΔCt methods. represents P <0.05.

4. Discussion

Osteoarthritis is the most common degenerative disease worldwide that affects small and large joints. Osteoarthritis often affects the whole joint tissue, including synovium, subchondral bone, and cartilage [12]. Synovitis is a common feature of osteoarthritis, even in early diseases. Synovial proliferation and tissue hypertrophy are significant in advanced osteoarthritis. There is increasing evidence that synovitis plays a key role in the pathogenesis and progression of osteoarthritis. Further understanding of the molecular and cellular variability of osteoarthritis-related synovitis can provide insight into the etiology of arthritis [4]. GWAS studies help to determine the important role of genetic factors in the risk of osteoarthritis. Understanding the key genes that influence the onset and progression of osteoarthritis will be necessary for the development of early treatment of the disease.

In our study, the analysis of gene expression profiling revealed the hub genes and pathways associated with osteoarthritis and enabled the identification of targets for therapeutic strategy. Bioinformatics methods are applied to analyze the raw data, and we identify 229 DEGs in osteoarthritis synovial membranes compared with normal synovial membranes, including 145 upregulated DEGs and 84 downregulated DEGs. Next, all DEGs were classified into three groups by GO terms using the online database and further clustered based on gene functions and signaling pathways, respectively. The PPI network of DEGs was constructed by STRING database, and the top ten hub genes were obtained: IL8, ADIPOQ, CXCR1, CXCL1, LIPE, FPR1, FABP4, SLC2A4, FASN, and IL6. The statistical result of the verification experiment showed that the expression level of ADIPOQ, IL6, and CXCR1 was significantly increased in knee osteoarthritis synovial membranes (p < 0.05). Studying these hub genes may contribute to the early diagnosis and treatment of osteoarthritis.

IL-6 plays a crucial role in chronic inflammation and the expression levels of IL-6 increase in human inflammatory diseases. In acute and chronic inflammation, IL-6 secreted into the serum and induced a transcriptional inflammatory response through interleukin 6 receptor alpha [13]. IL-6 expression was increased in osteoarthritis synovial membranes and synovial fluid, histone hyperacetylation, and DNA hypomethylation in the promoter of IL-6 gene were observed in osteoarthritis synovial fibroblasts compared with healthy synovial fibroblasts. Overexpression of IL-6 in synovial fibroblasts of osteoarthritis was suppressed through decrease histone acetylation and overexpression DNA methylation [14]. In our study, IL6 was significantly increased in knee osteoarthritis synovial membranes. The KEGG pathway enrichment showed that IL-6 was enriched in rheumatoid arthritis, cytokine-cytokine receptor interaction, and TNF signaling pathway [15, 16]. Previous studies showed that these pathways are involved in osteoarthritis synovial hyperplasia. Inflammatory cytokines, including IL-1, TNF, and IL-6, play a key role in osteoarthritis, and IL-6 is considered to be the key cytokine; its effect is mainly based on promoting the formation of osteoclasts and bone resorption while synergism with IL-1 and TNF [15]. In our previous research, the expression levels of IL6 and IL-1 were significantly increased in osteoarthritis synovial membranes [17]. Chemokine is a small secretory molecule, which mainly plays the role of chemotactic immune cells. The chemokine receptors CXCR1 and CXCR2 play a role in neutrophil-dependent injury and mediating neutrophil recruitment in inflammatory disease [18]. CXCR1/CXCR2 play a significant role in regulating tissue inflammation in a mouse model of osteoarthritis [19]. Cytokine networks in osteoarthritis may be the future research direction in therapeutic strategies. Different from microarray analysis results, ADIPOQ was significantly increased in knee osteoarthritis synovial membranes in our study. Adiponectin is a specific protein expressed in adipose tissue exclusively. The previous study found that SNP rs182052 is significantly associated with susceptibility to knee osteoarthritis in the Chinese population [20], and the interaction between rs1501299 (ADIPOQ) and rs662 (PON1) gene polymorphisms may play an important role in the development of osteoarthritis [21]. And the expression level of adiponectin was significantly higher in osteoarthritis patients than in controls [22]. In a word, ADIPOQ gene mutation may be associated with an increased risk of knee osteoarthritis.

In summary, using integrated bioinformatical analysis and qRT-PCR validation, we have identified the hub genes and related pathways, and these findings have the potential to be used as biomarkers and targets for osteoarthritis early diagnosis and treatment. There are still many limitations in our research: small sample size and lack of further experiments. To confirm our analysis results, more experimental studies with a larger sample are needed to confirm our study.

Data Availability

The microarray data used to support the findings of this study have been deposited in the GEO database (dataset ID:GSE89408), which have been cited.

Ethical Approval

Ethical number: this study was approved by the ethics committee of the Second Hospital of Jilin University (NO. 2018-082).

Conflicts of Interest

The authors have declared that no conflicts of interest exist.

Authors’ Contributions

Guizhen Zhang and Yang Song designed and supervised the experiments. Zhaoyan Li and Lei Zhong performed the experiments and analysis work. Zhenwu Du, Gaoyang Chen, and Jing Shang analyzed the data. Qingyu Wang and Qiwei Yang contributed to analysis tools.

Acknowledgments

This work was supported by the project supported by the National Natural Science Foundation of China (Grant No. 81702195), the Training Program of Outstanding Doctoral Student by Norman Bethune Health Science Center of Jilin University (No. 470110000646), the Graduate Innovation Fund of Jilin University (No. 101832018C075), the Science and Technology Department of Jilin Province of China (No. 20180520111JH), and the Second Hospital of Jilin University (No. KYPY2018-02).

References

M. Cross, E. Smith, D. Hoy et al., “The global burden of hip and knee osteoarthritis: estimates from the global burden of disease 2010 study,” Annals of the Rheumatic Diseases, vol. 73, no. 7, pp. 1323–1330, 2014.

K. E. Barbour, C. G. Helmick, M. Boring, and T. J. Brady, “Vital signs: prevalence of doctor-diagnosed arthritis and arthritis-attributable activity limitation – United States, 2013-2015,” Morbidity and Mortality Weekly Report (MMWR), vol. 66, no. 9, pp. 246–253, 2017.

E. J. Kubosch, G. Lang, D. Fürst et al., “The potential for synovium-derived stem cells in cartilage repair,” Current Stem Cell Research & Therapy, vol. 13, no. 3, pp. 174–184, 2018.

C. R. Scanzello and S. R. Goldring, “The role of synovitis in osteoarthritis pathogenesis,” Bone, vol. 51, no. 2, pp. 249–257, 2012.

L. J. Sandell, “Etiology of osteoarthritis: genetics and synovial joint development,” Nature Reviews Rheumatology, vol. 8, no. 2, pp. 77–89, 2012.

R. Edgar and A. Lash, The Gene Expression Omnibus (GEO): A Gene Expression and Hybridization Repository, National Center for Biotechnology Information, 2002.

Y. Guo, A. M. Walsh, U. Fearon et al., “CD40L-dependent pathway is active at various stages of rheumatoid arthritis disease progression,” The Journal of Immunology, vol. 198, no. 11, pp. 4490–4501, 2017.

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, p. e47, 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.

D. Szklarczyk, J. H. Morris, H. Cook et al., “The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible,” Nucleic Acids Research, vol. 45, no. 1, pp. D362–D368, 2017.

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, Supplement 4, no. S4, p. S11, 2014.

J. Martel-Pelletier, A. J. Barr, F. M. Cicuttini et al., “Osteoarthritis,” Nature Reviews Disease Primers, vol. 2, Article ID 16072, 2016.

C. Gabay, “Interleukin-6 and chronic inflammation,” Arthritis Research & Therapy, vol. 8, no. S2, p. S3, 2006.

F. Yang, S. Zhou, C. Wang et al., “Epigenetic modifications of interleukin-6 in synovial fibroblasts from osteoarthritis patients,” Scientific Reports, vol. 7, no. 1, article 43592, 2017.

P. Wojdasiewicz, Ł. A. Poniatowski, and D. Szukiewicz, “The role of inflammatory and anti-inflammatory cytokines in the pathogenesis of osteoarthritis,” Mediators of Inflammation, vol. 2014, Article ID 561459, 19 pages, 2014.

F. Xue, C. Zhang, Z. He, L. Ding, and H. Xiao, “Analysis of critical molecules and signaling pathways in osteoarthritis and rheumatoid arthritis,” Molecular Medicine Reports, vol. 7, no. 2, pp. 603–607, 2013.

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.

J. Sherwood, J. Bertrand, G. Nalesso et al., “A homeostatic function of CXCR2 signalling in articular cartilage,” Annals of the Rheumatic Diseases, vol. 74, no. 12, pp. 2207–2215, 2015.

F. M. Coelho, V. Pinho, F. A. Amaral et al., “The chemokine receptors CXCR1/CXCR2 modulate antigen-induced arthritis by regulating adhesion of neutrophils to the synovial microvasculature,” Arthritis & Rheumatism, vol. 58, no. 8, pp. 2329–2337, 2014.

L. Jiang, X. Zhu, J. Rong et al., “Obesity, osteoarthritis and genetic risk: the rs182052 polymorphism in the ADIPOQ gene is potentially associated with risk of knee osteoarthritis,” Bone & Joint Research, vol. 7, no. 7, pp. 494–500, 2018.

J. Fernández-Torres, G. A. Martínez-Nava, Y. Zamudio-Cuevas, K. Martínez-Flores, and R. Espinosa-Morales, “Epistasis between ADIPOQ rs1501299 and PON1 rs662 polymorphisms is potentially associated with the development of knee osteoarthritis,” Molecular Biology Reports, vol. 46, no. 2, pp. 2049–2058, 2019.

Q. Tang, Z. C. Hu, L. Y. Shen, P. Shang, H. Z. Xu, and H. X. Liu, “Association of osteoarthritis and circulating adiponectin levels: a systematic review and meta-analysis,” Lipids in Health and Disease, vol. 17, no. 1, article 189, 2018.