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

Li, Zhaoyan; Zhong, Lei; Du, Zhenwu; Chen, Gaoyang; Shang, Jing; Yang, Qiwei; Wang, Qingyu; Song, Yang; Zhang, Guizhen

doi:https://doi.org/10.1155/2019/8340573

BioMed Research International

On this page

Abstract Introduction Materials and Methods Results Discussion Data Availability Ethical Approval Conflicts of Interest Authors’ Contributions Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2019 | Article ID 8340573 | https://doi.org/10.1155/2019/8340573

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

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

Academic Editor: Magali Cucchiarini

Received11 Mar 2019

Revised19 Mar 2019

Accepted08 Apr 2019

Published16 Apr 2019

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.

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.

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.

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

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.

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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
C. R. Scanzello and S. R. Goldring, “The role of synovitis in osteoarthritis pathogenesis,” Bone, vol. 51, no. 2, pp. 249–257, 2012.
View at: Publisher Site | Google Scholar
L. J. Sandell, “Etiology of osteoarthritis: genetics and synovial joint development,” Nature Reviews Rheumatology, vol. 8, no. 2, pp. 77–89, 2012.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
J. Martel-Pelletier, A. J. Barr, F. M. Cicuttini et al., “Osteoarthritis,” Nature Reviews Disease Primers, vol. 2, Article ID 16072, 2016.
View at: Publisher Site | Google Scholar
C. Gabay, “Interleukin-6 and chronic inflammation,” Arthritis Research & Therapy, vol. 8, no. S2, p. S3, 2006.
View at: Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar

Copyright

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.

PDF Download Citation

Download other formats

Order printed copies

Views

2644

Downloads

1540

Citations