Identification of Key Genes and Pathways Associated with Differences of Subchondral Bone in Osteoarthritis

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

Download PDF

Research article

Identification of Key Genes and Pathways Associated with Differences of Subchondral Bone in Osteoarthritis

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

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

Osteoarthritis (OA) is a common chronic disease worldwide. Subchondral bone is an important pathological change in OA and responds more rapidly to adverse loading and events compared to cartilage. However, the pathogenic genes and pathways of subchondral bone are largely unclear.

Objective

This study aimed to identify signature differences in genes involved in knee lateral tibial (LT) and medial tibial (MT) plateaus of subchondral bone tissue while exploring their potential molecular mechanisms via bioinformatics analysis.

Methods

First, the gene expression data of GSE51588 was downloaded from the GEO database. Differentially expressed genes (DEGs) between knee LT and MT were identified, and functional enrichment analyses were performed. Then, a protein-protein interactive network was constructed in order to acquire the hub genes, and modules analysis was conducted using STRING and Cytoscape for further analysis. The enriched hub genes were queried in DGIdb database to find suitable drug candidates in OA.

Results

A total of 202 DEGs (112 upregulated genes and 84 downregulated genes) were determined. In the PPI network, ten hub genes were identified. Five significant modules were identified using the MCODE plugin unit. Functional enrichment analysis revealed the most important signaling pathways. Six of the ten hub genes were targetable by a total of 35 drugs, suggesting their possible therapeutic use for OA .

Conclusions

The identified hub genes and functional enrichment pathways were implicated in the development and progression of subchondral bone in OA, thus improving our understanding of OA and offering molecular targets for future therapeutic modalities.

Orthopedics

Bioinformatics analysis

osteoarthritis

subchondral bone

differentially expressed gene

protein-protein interaction network

microarray

Osteoarthritis (OA) is one of the most common forms of arthritis in the elderly population, which affects more than 6‑32% of people worldwide [1]. In previous studies, researchers have primarily focused on cartilage tissue and chondrocytes in regard to OA treatment and their associated mechanisms [2, 3]. Due to rapid advances in medical technologies, synovial tissue as well as subchondral bone has been further investigated[4]. Unfortunately, existing therapies are unable to prevent or reverse the development of OA apart from providing pain relief [5]. Eventually, patients with OA in advanced stages commonly receive joint replacement surgery. Hence, there is an urgent need to explore potential biomarkers and therapeutic targets for OA.

Subchondral bone located under articular cartilage is an important part of the joint as it is a composite tissue composed of mineral, organic matrix (non-collagenous proteins, collagen, and lipids), cells, and water [6]. Anatomically, subchondral bone consists of both the underlying subchondral trabecular bone and a layer of corticalized subchondral bone plate. Due to subchondral bone’s specialized structure, it may effectively spare the overlying articular cartilage by dispersing axial loads across the joint [7]. The major function of subchondral bone is to absorb stress and maintain the shape of the joint [8]. According to Goldring (2010), subchondral bone responds more rapidly to adverse loading and events compared to cartilage [9]. By visualizing the whole knee’s structure via magnetic resonance imaging (MRI), the severity of knee OA is associated with increased quantities of tibial subchondral bone [10]. As subchondral sclerosis increases, a reduction in the density of the underlying subchondral bone and flexibility of the osteochondral junction is present. The texture of subchondral bone may generally be regarded as an important biomarker to predict the progression of OA [11]. Subchondral bone is currently considered a challenge in view of achieving osteochondral regeneration and self-repair [12]. The aforementioned studies provide an important foundation for further exploration of the development of OA. However, the mechanism of aberrant subchondral bone formation and articular cartilage degeneration remains unclear in OA.

Due to strides in biomedicine, numerous studies have shown that the occurrence and development of OA may be mediated by multiple genes and signaling pathways [13]. In addition, microarray technologies have been used to perform multiple gene expression profiling studies on human knee OA tissue. Integrated bioinformatics techniques combined with expression profiling may offer additional advantages compared to independent studies [14]. The majority of biomedical studies using gene expression microarray profiling thus far have focused on human articular cartilage [15], synovium[16], or synovial fibroblasts [17]. Consequently, none have been performed on human subchondral bone tissue, particularly in regard to LT and MT.

This study attempts to identify the biomarkers, pathways and therapeutic targets involved in LT and MT of knee subchondral bone tissue via bioinformatics analysis and microarray gene expression profiling in order to elucidate the subchondral bone tissue mechanism involved in OA and develop more effective treatment options.

2.1. Microarray Data Information.

The NCBI Gene Expression Omnibus (GEO) is a public genomics data repository (http://www.ncbi.nlm.nih.gov/geo) that stores gene expression datasets, platform records and original series [18]. Gene expression profile dataset GSE51588 [19] was obtained from the GEO database, which was expressed on the GPL13497 platform [Agilent-026652 Whole Human Genome Microarray 4 × 44K v2 (Probe Name version), Agilent Technologies Manufacturer]. The GSE51588 dataset, which was the overlying articular cartilage or synovium cells, contained total RNA from the human knee (40 OA and 10 non-OA) from regions of LT and MT. Only 40 samples with OA were selected in our study for further investigation.

2.2. Data Processing and Identification of DEGs.

The common statistical method entailed in DEG analysis is to select genes at significantly different expression levels between samples [20]. We used the R software (version 3.6.1) and limma packages [21] to identify DEGs between knee lateral tibial (LT) and medial tibial (MT) plateaus in patients with OA. Then, the P-values were adjusted for comparisons using the false discovery rate (FDR) of the Benjamini and Hochberg (BH) test [22] in the limma package. The criteria for DEG was |log2-fold change (FC)| ≥2 and adjusted P-values < 0.01.

2.3. Functional and pathway enrichment analysis

The Database for Annotation, Visualization, and Integrated Discovery (DAVID; https://david.ncifcrf.gov/) [23] integrates biological data and analysis tools that provide systematic information for the list of DEGs. Gene Ontology (GO, www.geneontology.org) is widely used in the annotation of cell components (CCs), molecular functions (MFs), and Biological processes (BPs) of DEGs. The KEGG (Kyoto Encyclopedia of Genes and Genomes; www.genome.ad.jp/kegg) is a database resource used to extract pathway information from molecular interaction networks [24]. In the present study, GO and KEGG pathway enrichment analysis of DEGs were conducted using the DAVID online tool. The p-value <0.05 and gene count >2 were set as the cutoff point.

2.4. Construction of the PPI network analysis

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (http://www.string-db.org/) is an online biological database and website designed to construct the PPI networks in molecular biology [25]. The DEGs (species: Homo) were mapped to STRING(PPI score >0.4), and then Cytoscape software 3.7.1 (http://www.cytoscape.org/) was used to visualize the PPI network. Afterward, the hub genes were identified according to degree of connectivity.

2.5. Functional module analysis of PPI network

Functional module analysis of PPI networks was performed using the molecular complex detection (MCODE) to cluster the algorithm accompanied with Cytoscape [26]. The analysis parameters were set by default. Additionally, the module analysis results were executed for the functional and pathway enrichment analyses in each module using DAVID, with p<0.05 set as the threshold indicating a statistically significant difference.

2.6. Drug‑gene interactions.

This study employed the Drug-gene Interaction Database (DGIdb; http://www.dgidb.org), which was used as a potential target to search for existing drugs [27] in order to drug‑gene interactions between the final hub genes. Such drugs targeting genes or pathways may reveal potential treatments of OA.

3.1. Identification of DEGs

In the present study, DEGs were identified from 20 paired overlying articular cartilage tissues via Limma. By comparing LT and MT in patients suffering from OA, a total of 202 DEGs were identified, including 118 upregulated and 84 downregulated genes (Table 1). The red represents upregulated genes while blue represents downregulated genes in the volcano plot (Fig. 1). The heat map regarding the top 30 upregulated and downregulated DEGs are shown in Fig. 2.

Table 1

Two hundred and two EDGs were identified in GSE51588, including 118 Up-regulated Genes and 84 Down-regulated Genes in the patients with OA.a
DEGs	Gene name
Up-regulated genes	LEP, APOB, ADH1B, CCL8, TUSC5, KLB, SLC19A3, NPR3, PLIN1, AQP7P1, LOC283392, PCK1, MYO16, PLIN4, LIPE, EYS, ADH1A, SHISA6, SORBS1, FOSB, DGAT2, PKP2, MGST1, G0S2, NPY1R, COPG2IT1, AQP7P3, FMO2, AQP7, CA4, ITIH5, CYP4B1, ANKRD30A, PLEKHA6, MRAP, ADIPOQ, CIDEC, OLFM4, CHRDL1, TMEM132C, ACADL, KIAA0125, CIDEA, ADRB1, NTRK2, PLXNA4, FABP4, ADH1C, MLXIPL, FOS, CXCL2, CIDECP, GPD1, CYP4 × 1, EGFEM1P, LMO3, RGS6, VWF, SLC2A4, DMRT2, MYZAP, LPL, CPA1, CCL3, ACVR1C, P2RY14, HCAR3, ANKRD20A5P, ANKRD20A2, ANKRD20A9P, FAM70A, EPB41L4B, ATP1A2, AQPEP, CEBPA, MAOB, FREM1, LOC401847, PNPLA2, TIMP4, MAPT, CHRDL2, PDE3B, MAOA, IRF4, SPOCK3, GPAM, CADM2, LEPR, C14orf91, 1-Mar, PRKAR2B, PPP1R14A, RDH5, KCTD8, PLIN5, AMPD1, TMC5, TTPA, TNFRSF13B, LOC152286, TMEM56, CCL4, MESP1, RSPO3, IGLL5, SPON1, TNS1, ANKRD20A8P, ELMOD3, PTGER3, C14orf180, MYOM1, ANTXRL, SLC24A3, TSPAN8, AOC3, PDE8B
Down-regulated genes	STMN2, IL11, POSTN, APCDD1L, GPR179, HMBOX1, LOC389043, LRRC2, ASPN, COL6A3, FAM159B, CD86, OR52E8, PPP6R1, TUBB3, COL3A1, MAGEB6, KRT79, C21orf37, CA12, BICC1, RAP1GAP, OGN, ATXN3L, CCL18, DIO2, PPEF1, ST6GAL2, GDF6, TNFSF11, LOC100287221, PBK, DSG3, PROM2, WNT16, RPA4, CRLF1, THBS4, CTHRC1, KCNK2, SBSN, MGC4294, SOX11, C13orf30, GRP, XAGE-4, PAX1, CENPA, VPS18, FNDC1, ANLN, SLC13A5, C6orf138, MYO3A, CEP55, LOC283731, MMP13, LHX2, ELMOD1, COL17A1, THY1, RALYL, SLC16A10, VSNL1, FAP, CXADR, PLEK2, BIRC5, JAKMIP2, DNER, SPEN, PTGES, PRSS35, RGS13, C9orf53, GNG4, ECM1, TNNT3, ALPK2, WDR33, C1QTNF6, CPZ, OLFML2B, MGC39584
DEG, differentially expressed gene ^aThe differentially expressed genes are listed from the smallest to the largest adjusted p value.

3.2. Functional and pathway enrichment analysis

To further explore the potential targets of these DEGs between LT and MT, we performed GO and KEGG pathway analyses (Figs. 3 and 4). As illustrated in Fig. 3, we classified GO analysis into 3 functional categories: including top ten of BP (cellular response to tumor necrosis factor, response to tumor necrosis factor, temperature homeostasis, etc.), top three of CC (collagen − containing extracellular matrix, lipid droplet, collagen trimer), and the top ten of MF (receptor-ligand activity, extracellular matrix structural constituent, cytokine activity, etc.). The top ten of KEGG pathway (PPAR signaling pathway, Tyrosine metabolism, and Cytokine-cytokine receptor interaction, etc.) were described in Fig. 4.

We also respectively show annotation of the up-regulated genes as well as down-regulated genes (Tables 2 and 3). In the BP group, the up-regulated DEGs were mainly enriched for genes involved in the triglyceride metabolic process, chemical homeostasis, and acylglycerol metabolic process. The downregulated DEGs were enriched genes that participated in epithelial cell proliferation, negative regulation of developmental processes, and locomotion. In the CC group, the upregulated DEGs were mainly enriched for genes associated with lipid particles, intrinsic components of the plasma membrane, and integral components of the plasma membrane. Moreover, the downregulated DEGs were enriched for genes associated with the proteinaceous extracellular matrix, extracellular matrix, and collagen trimer. In the MF group, the up-regulated DEGs were mainly enriched from genes associated with protein homodimerization activity, protein dimerization activity, and alcohol dehydrogenase activity; and the down-regulated DEGs were enriched for genes encoding proteins for receptor binding, integrin binding, and cytokine receptor binding. In the KEGG pathway group (P < 0.005), the up-regulated DEGs were enriched for genes associated with a total of nine terms (Regulation of lipolysis in adipocytes, PPAR signaling pathway, Tyrosine metabolism, etc.) but down-regulated DEGs only one term (protein digestion and absorption). (Table 4)

Table 2

The Top Three Pathways in GO Enrichment Analysis of Up-regulated DEGs
Category	Term	Description	P-Value	Count	Gene Set
BP	GO:0006641	triglyceride metabolic process	6.40E-09	10	GPD1, LPL, PLIN5, APOB, DGAT2, FABP4, PNPLA2, GPAM, LIPE, PCK1
BP	GO:0048878	chemical homeostasis	1.27E-08	24	CEBPA, LPL, CCL3, PTGER3, LEPR, MLXIPL, TTPA, CCL8, AQP7, PNPLA2, ATP1A2, ACADL, ADIPOQ, PCK1, LEP, PRKAR2B, APOB, SLC2A4, DGAT2, SLC24A3, PKP2, FABP4, GPAM, AQP7P3
BP	GO:0006639	acylglycerol metabolic process	2.13E-08	10	GPD1, LPL, PLIN5, APOB, DGAT2, FABP4, PNPLA2, GPAM, LIPE, PCK1
CC	GO:0005811	lipid particle	2.48E-10	10	PLIN5, DGAT2, PLIN1, PLIN4, CIDEA, FABP4, G0S2, PNPLA2, LIPE, CIDEC
CC	GO:0031226	intrinsic component of plasma membrane	0.001396	22	PLXNA4, PTGER3, KLB, CADM2, TNFRSF13B, AQP7, TSPAN8, SLC19A3, NPY1R, ATP1A2, NPR3, ACVR1C, SHISA6, ADRB1, SLC2A4, SORBS1, SLC24A3, P2RY14, NTRK2, CA4, HCAR3, AQP7P3
CC	GO:0005887	integral component of plasma membrane	0.002076	21	PLXNA4, PTGER3, KLB, CADM2, TNFRSF13B, AQP7, TSPAN8, SLC19A3, NPY1R, ATP1A2, NPR3, ACVR1C, SHISA6, ADRB1, SLC2A4, SORBS1, SLC24A3, P2RY14, NTRK2, HCAR3, AQP7P3
MF	GO:0042803	protein homodimerization activity	1.22E-05	16	CEBPA, GPD1, OLFM4, CADM2, MAOB, MLXIPL, DMRT2, CIDEA, NPR3, ADIPOQ, VWF, DGAT2, NTRK2, MYOM1, MGST1, AOC3
MF	GO:0046983	protein dimerization activity	6.55E-05	19	CEBPA, GPD1, OLFM4, CADM2, MAOB, MLXIPL, DMRT2, CIDEA, NPR3, ADIPOQ, FOS, VWF, ADRB1, DGAT2, NTRK2, MYOM1, MESP1, MGST1, AOC3
MF	GO:0004024	alcohol dehydrogenase activity	4.62E-04	3	ADH1C, ADH1B, ADH1A
GO: gene ontology; DEGs: differentially expressed genes

Table 3

The Top Three Pathways in GO Enrichment Analysis of Down-regulated DEGs
Category	Term	Description	P-Value	Count	Gene Set
BP	GO:0050673	epithelial cell proliferation	3.26E-04	8	WNT16, TNFSF11, RAP1GAP, SOX11, FAP, BIRC5, ECM1, THBS4
BP	GO:0051093	negative regulation of developmental process	4.23E-04	11	ASPN, CD86, RAP1GAP, STMN2, SOX11, COL3A1, POSTN, ECM1, KCNK2, THBS4, THY1
BP	GO:0040011	locomotion	0.00146	15	CTHRC1, COL3A1, ANLN, POSTN, CXADR, ECM1, CCL18, THY1, OGN, TNFSF11, DNER, LHX2, FAP, TUBB3, THBS4
CC	GO:0005578	proteinaceous extracellular matrix	2.21E-08	13	ASPN, CTHRC1, WNT16, OLFML2B, COL3A1, POSTN, ECM1, MMP13, CPZ, OGN, COL17A1, COL6A3, THBS4
CC	GO:0031012	extracellular matrix	1.88E-07	14	ASPN, CTHRC1, WNT16, OLFML2B, COL3A1, POSTN, ECM1, MMP13, CPZ, OGN, COL17A1, SBSN, COL6A3, THBS4
CC	GO:0005581	collagen trimer	3.81E-06	7	CTHRC1, COL17A1, C1QTNF6, COL3A1, COL6A3, MMP13, WDR33
MF	GO:0005102	receptor binding	1.86E-04	17	CTHRC1, WNT16, GDF6, COL3A1, CRLF1, CXADR, ECM1, CCL18, THY1, IL11, GRP, OGN, CD86, DNER, FAP, THBS4
MF	GO:0005178	integrin binding	9.08E-04	5	FAP, COL3A1, CXADR, THBS4, THY1
MF	GO:0005126	cytokine receptor binding	0.005027	6	TNFSF11, GDF6, CRLF1, ECM1, CCL18, IL11
GO: gene ontology; DEGs: differentially expressed genes

Table 4

The Top Pathways in KEGG Enrichment Analysis of DEGs (p < 0.005)
Category	Term	Description	P-Value	Count	Gene Set
Up-regulated DEGs	hsa04923	Regulation of lipolysis in adipocytes	8.96E-09	9	PTGER3, ADRB1, PLIN1, PDE3B, FABP4, AQP7, PNPLA2, NPY1R, LIPE
	hsa03320	PPAR signaling pathway	7.45E-07	8	LPL, SORBS1, PLIN1, FABP4, AQP7, ACADL, ADIPOQ, PCK1
	hsa00350	Tyrosine metabolism	6.70E-06	6	MAOA, MAOB, ADH1C, ADH1B, ADH1A, AOC3
	hsa00982	Drug metabolism-cytochrome P450	1.33E-05	7	FMO2, MAOA, MAOB, ADH1C, ADH1B, ADH1A, MGST1
	hsa04024	cAMP signaling pathway	8.28E-04	8	FOS, PTGER3, ADRB1, PDE3B, ATP1A2, NPY1R, HCAR3, LIPE
	hsa04920	Adipocytokine signaling pathway	0.00212	5	LEP, SLC2A4, LEPR, ADIPOQ, PCK1
	hsa04152	AMPK signaling pathway	0.002646	6	LEP, SLC2A4, LEPR, ADIPOQ, LIPE, PCK1
	hsa00071	Fatty acid degradation	0.004229	4	ADH1C, ADH1B, ADH1A, ACADL
	hsa04910	Insulin signaling pathway	0.004343	6	PRKAR2B, SLC2A4, SORBS1, PDE3B, LIPE, PCK1
Down-regulated DEGs	hsa04974	Protein digestion and absorption	0.003796	4	COL17A1, SLC16A10, COL3A1, COL6A3
KEGG, Kyoto Encyclopedia of Genes and Genomes

3.3. Construction of the PPI network analysis

Based on the data in the STRING database, we constructed a PPI network through the Cytoscape software, which contained 130 nodes and 362 edges (Fig. 5). Among the 130 genes, the top 10 hub genes were identified according to degree connectivity, including LEP, PNPLA2, LIPE, PLIN1, ADIPOQ, FABP4, LPL, SLC2A4, DGAT2, and CIDEC. LEP demonstrated the highest degree score (score = 31). (Table 5)

Table 5

The top 10 of DEGs with higher degree centrality
Node	DC	DEGs	BC	CC	NC
LEP	31	Up-regulated	0.369552	0.420863	12
PNPLA2	24	Up-regulated	0.058373	0.336207	12.91667
LIPE	24	Up-regulated	0.03605	0.346154	13.45833
PLIN1	21	Up-regulated	0.035155	0.345133	14.28571
ADIPOQ	20	Up-regulated	0.057178	0.354545	14.95
FABP4	20	Up-regulated	0.030822	0.345133	15
LPL	18	Up-regulated	0.015573	0.336207	16.5
SLC2A4	18	Up-regulated	0.032263	0.343109	15.27778
DGAT2	18	Up-regulated	0.014033	0.334286	15.94444
CIDEC	16	Up-regulated	0.010125	0.323204	15.625
DC: degree centrality, DEGs: differentially expressed genes, BC: betweenness centrality, CC: closeness centrality, NC: Neighborhood connectivity

3.4. Functional module analysis of PPI networks

To further analyze the interactions of proteins, 5 modules were detected using the Cytoscape plugin MCODE (Table 6 and Fig. 5). In addition, functional enrichment analyses for these modules were performed. The pathway enrichment analysis showed that Module 1 was mainly associated with PPAR signaling pathway, AMPK signaling pathway, and Regulation of lipolysis in adipocytes. Module 2 was found to be mainly associated with Chemokine signaling pathway, cAMP signaling pathway, and Neuroactive ligand-receptor interaction. Furthermore, Module 3 was mainly associated with Tyrosine metabolism, Drug metabolism - cytochrome P450, and Fatty acid degradation.

Table 6.
Basic Characteristics of Top Five Module
Module	Score	Node	Edges	Gene names
Cluster 1	13	14	85	APOB, MLXIPL, PNPLA2, LPL, PLIN1, LIPE, CIDEA, ADIPOQ, FABP4, SLC2A4, DGAT2, LEP, CIDEC, PCK1
Cluster 2	7	7	21	HCAR3, PTGER3, CCL4, GNG4, P2RY14, NPY1R, CXCL2
Cluster 3	5	5	10	MAOA, ADH1A, MAOB, ADH1B, ADH1C
Cluster 4	5	5	10	COL3A1, OGN, POSTN, COL6A3, ASPN
Cluster 5	5	5	10	BIRC5, CEP55, ANLN, PBK, CENPA
Parameter setting: Find Clusters: In Whole Network; Degree Cutoff: 2; Cluster Finding: Haircut; Node Score Cutoff: 0.2; K-core: 2; Max. Depth: 100

Table 7

Candidate drugs targeting genes with OA.
Number	Drug	Gene	Approved by FDA	Approved use in OA	Drug categories	Approved use
1	METHADONE	LEP	Yes	No	Synthetic opioid analgesic	Management of pain
2	GEMFIBROZIL	LEP	Yes	No	Fibric acid agent	Hyperlipidemia, CHD
3	FOLIC ACID	LEP	Yes	No	Folate or Vitamin B9	Folic acid deficiency
4	BROMOCRIPTINE	LEP	Yes	No	Stimulating dopaminergic receptors	PD,hyperprolactinemia
5	LOVASTATIN	LEP	Yes	No	Statin medications	Cardiovascular disease,Dyslipidemia
6	LAMIVUDINE	LEP	Yes	No	Zalcitabine analog	HIV, HBV Infections
7	MIRTAZAPINE	LEP	Yes	No	Piperazino-azepine agent	Major depressive disorder
8	METHIMAZOLE	LEP	Yes	No	Hormone Antagonists	Thyroid disease
9	FELODIPINE	LEP	Yes	No	Long-acting dihydropyridine CCB	Hypertension
10	ISOPROTERENOL	LIPE	Yes	No	Isopropyl analog	Bronchodilator and heart stimulant
11	FENOFIBRATE	ADIPOQ	Yes	No	Fibric acid derivative	Dyslipidemia
12	NEVIRAPINE	ADIPOQ	Yes	No	NNRTI	Ongoing HIV-1
13	SPIRONOLACTONE	ADIPOQ	Yes	No	Potassium sparing diuretic	III-IV heart failure, Hepatocirrhosis
14	DIAZOXIDE	LPL	Yes	No	A benzothiadiazine derivative	hypertensive emergencies
15	MIFEPRISTONE	LPL	Yes	No	Glucocorticoid receptor antagonist	Endocrine diseases
16	LABETALOL	LPL	Yes	No	Adrenergic Antagonists	Hypertension
17	TRIAMCINOLONE	LPL	Yes	Yes	Adrenal Cortex Hormones	Various inflammatory conditions
18	ISOTRETINOIN	LPL	Yes	No	Retinoid derivative	Recalcitrant nodular acne
19	ACARBOSE	LPL	Yes	No	A glucosidase inhibitor	Diabetes
20	ASPARAGINASE	LPL	Yes	No	Amidohydrolases	ALL
21	NADOLOL	LPL	Yes	No	Adrenergic Antagonists	Angina pectoris, hypertension
22	NYSTATIN	SLC2A4	Yes	No	Polyene antifungal drug	Fungal infections
23	ETOPOSIDE	SLC2A4	Yes	No	Antineoplastic Agents	Malignancies
24	PENICILLAMINE	SLC2A4	Yes	No	Immunosuppressive Agents	Wilson's disease
25	GLYBURIDE	SLC2A4	Yes	No	Generation sulfonylurea	Diabetes mellitus
26	INDINAVIR	SLC2A4	Yes	No	Antiviral Agents	HIV infection
27	DIPYRIDAMOLE	SLC2A4	Yes	No	Phosphodiesterase inhibitor	Thromboembolic complications
28	NEOMYCIN	SLC2A4	Yes	No	Aminoglycoside antibiotic drug	Bacterial infection
29	INDOMETHACIN	SLC2A4	Yes	Yes	NSAID	Rheumatoid arthritis
30	IMATINIB	SLC2A4	Yes	No	Antineoplastic Agents	Leukaemia
31	VERAPAMIL	SLC2A4	Yes	No	Phenylalkylamine CCB	Hypertension, cardiopathy
32	IRBESARTAN	SLC2A4	Yes	No	ARB	Hypertension, diabetic nephropathy
33	STREPTOZOTOCIN	SLC2A4	Yes	No	Antibiotics, Antineoplastic	Malignant neoplasms
34	COLCHICINE	SLC2A4	Yes	No	Alkaloid drug	Arthrolithiasis
35	NIACIN	DGAT2	Yes	No	Vitamin B Complex	Vitamin deficiencies
CHD: Coronary Heart Disease; PD: Parkinson's Disease; HIV: Human Immunodeficiency Virus; HBV: Hepatitis B Virus; CCB: calcium channel blocker; NNRTI: non-nucleoside reverse transcriptase inhibitor; ALL: Acute Lymphoblastic Leukaemias; NSAID: non-steroidal anti-inflammatory drug; ARB: angiotensin receptor blocker

3.5 Drug‑gene interactions.

Using the10 hub genes the PPI network analysis revealed the potential targets. Finally, only five genes were explored as potential drugs: LEP (9 drugs), ADIPOQ (3 drugs), LPL (8 drugs), SLC2A4 (13 drugs), and LIPE and DGAT2 (1 drug each). A list of the 35 drugs was compiled, which were all approved for use in humans by the FDA, satisfied the standards for being possible drug treatments for OA (Table III and Fig. 6).

Despite demonstrated risk factors such as obesity, mechanical abnormality, genetic predisposition, and age, the precise pathogenesis underlying OA remains unclear; which leads to the present therapeutic program for OA is only primarily relieving symptoms not specific treatment. Several prior OA studies emphasized only synovium [28], articular cartilage [29], and meniscus [30] while ignoring the role of subchondral bone in OA development [31]. Recently, increasing evidence has indicated that the degeneration of articular cartilage was associated with abnormal metabolism of subchondral bone[32, 33]. Microarray technology and bioinformatics techniques have been widely used in predicting the potential therapeutic targets. This study attempted to further understand the molecular mechanisms involved in knee subchondral bone tissue of LT and MT using Microarray technology and bioinformatics.

To the best of our knowledge, this is the first study that identifies the key genes and pathways between LT and MT from knee subchondral bone tissue using bioinformatics analysis. Data was extracted from GSE51588, including 40 OA (20 LT and 20 MT from human knee subchondral bone tissue) and 10 non-OA samples (5 LT and 5 MT from human knee subchondral bone tissue). Only 40 OA samples were enrolled in the present study, and a total of 202 DEGs were identified between LT and MT of knee subchondral bone, including 119 upregulated and 84 downregulated genes. Moreover, 10 genes were further screened with higher degree scores (> 15) and the most significant expression changes, including LEP, PNPLA2, LIPE, PLIN1, ADIPOQ, FABP4, LPL, SLC2A4, DGAT2, and CIDEC. As numerous studies have demonstrated that co-expressed genes are frequently related to those with similar expression profiles that are also commonly involved in similar biological processes [34], to better reveal the functions of the identified DEGs, we performed GO and KEGG enrichment analyses as well as constructed of the PPI network analysis.

The KEGG enrichment analysis of the DEGs demonstrated that the hub genes PLIN1, FABP4, PNPLA2, and LIPE were mainly enriched in regulation of lipolysis in the adipocytes; LPL, PLIN1, FABP4, and ADIPOQ were mainly involved in PPAR signaling pathway, while LEP, SLC2A4, ADIPOQ, and LIPE were mainly enriched in AMPK signaling pathway. Adipocytes are derived from multipotent mesenchymal stem cells (MSCs), and various molecular events commit MSCs to the adipocyte lineage [35]. However, MSCs are treated as cellular sources for multiple tissue regeneration, which can differentiate into different cell lineages, like muscle, bone, cartilage, or fat [36]. These studies tested the importance in the pathway of regulation of lipolysis in adipocytes for OA. Furthermore, in a study on the overexpression of miR-10a-5p inducing human primary chondrocyte, they also enriched important terms including PPAR signaling pathway [37]. This finding is also consistent with the teams of Li [38] and Yi [39] those who performed enrichment analysis using the tissue of knee OA synovial membranes. The α subunit of AMPK exists in two isoforms in the mammalian, including AMPKα1 and AMPKα2, where both serve as catalytic subunits [40]. Increased NR4A1 acts on the AMPK/dynaminrelated protein 1 (Drp1) pathway to trigger mitochondrial fission, leading to mitochondrial permeability transition pore related chondrocyte cell death [41]. Such results are consistent with that of the present study and highlight the credibility of the enrichment analysis.

Moreover, PPI networks involved with DEGs were constructed, where the top 10 hub genes were observed to be LEP, PNPLA2, LIPE, PLIN1, ADIPOQ, FABP4, LPL, SLC2A4, DGAT2, and CIDEC. LEP has mitogenic activity on vascular endothelial cells and promotes angiogenesis [42]. In the early stage of OA, the high expression level of LEP may promote the phosphorylation of AKT [43], increasing bone remodeling via triggering the downstream signaling pathway. Prior studies have shown that skeletal muscle dysfunction and decline in mobility are associated with the PNPLA2/ATGL pathway [44]. LIPE is in white adipose tissue of retinyl ester hydrolase; LIPE-null mice have shown decreased levels of retinaldehyde, retinol and all-trans RA as well as increased levels of retinyl esters [45]. PLIN proteins in mammalian involved a conserved PAT domain (Plin1, Plin2 and Plin3) which appears to be important in lipid droplet targeting [46]. Early clinical studies have shown that the expression level of ADIPOQ are correlated with disease severity of knee OA [47]. ADIPOQ may be treated as a potential biomarker for knee damage in OA [48]. Six hub genes (LEP, LIPE, LPL, SLC2A4, FABP4,and ADIPOQ) corroborate the findings of Gu et al. (2019) who used weighted gene co-expression network analysis to identify the biomarkers of subchondral bone [49].

In the present study, a bioinformatics analysis was performed to investigate knee LT and MT plateaus of subchondral bone tissue in OA. The identified functional enrichment pathways and 10 hub genes were identified providing further insight into the underlying pathogenesis and facilitate the diagnosis and treatment of OA. However, a limitation was that the screened genes and pathways lacked additional experimentation. Therefore, further studies are required to confirm the corresponding results.

OA, Osteoarthritis; LT, lateral tibial;MT, medial tibial; GEO, Gene Expression Omnibus; DEGs, Differentially expressed genes; PPI, protein-protein interactive; DGIdb, Drug-Gene Interaction database; FDR, false discovery rate; BH, Benjamini and Hochberg; DAVID, Database for Annotation, Visualization, and Integrated Discovery; CCs, cell components; MFs, molecular functions; BPs, Biological processes; STRING, Search Tool for the Retrieval of Interacting Genes/Proteins; KEGG, Kyoto Encyclopedia of Genes and Genomes; MCODE, molecular complex detection

Authors’ contributions

RGY performed comparative analysis using bioinformatics tools. JYZ, EYF, HYW, YGZ, WLW, ZLL, LQL,YTH, and CLW participated in data analysis and discussion, RGY interpreted data and wrote the manuscript. YYZ organized and supervised the project. All authors read and approved the final manuscript. RGY was responsible for manuscript review.

Acknowledgments

Not applicable.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Availability of data and materials

The datasets generated during and/or analysed during the current study are available in the [Pubmed] repository, [doi:10.1186/ar4380]

Funding Information

No funding was received for this article.

McDougall C, Hurd K, Barnabe C (2017) Systematic review of rheumatic disease epidemiology in the indigenous populations of Canada, the United States, Australia, and New Zealand. Seminars in arthritis and rheumatism 46 (5):675-686. doi:10.1016/j.semarthrit.2016.10.010
Hashimoto S, Ochs RL, Komiya S, Lotz M (1998) Linkage of chondrocyte apoptosis and cartilage degradation in human osteoarthritis. Arthritis and rheumatism 41 (9):1632-1638. doi:10.1002/1529-0131(199809)41:9<1632::Aid-art14>3.0.Co;2-a
Schroeppel JP, Crist JD, Anderson HC, Wang J (2011) Molecular regulation of articular chondrocyte function and its significance in osteoarthritis. Histology and histopathology 26 (3):377-394. doi:10.14670/hh-26.377
Wang X, Hunter DJ, Jin X, Ding C (2018) The importance of synovial inflammation in osteoarthritis: current evidence from imaging assessments and clinical trials. Osteoarthritis and cartilage 26 (2):165-174. doi:10.1016/j.joca.2017.11.015
Seed SM, Dunican KC, Lynch AM (2011) Treatment options for osteoarthritis: considerations for older adults. Hospital practice (1995) 39 (1):62-73. doi:10.3810/hp.2011.02.375
Boskey AL (2010) The Biochemistry of Bone: Composition and Organization. In: Osteoporosis in Men. Elsevier, pp 3-13
Brama PA, Holopainen J, van Weeren PR, Firth EC, Helminen HJ, Hyttinen MM (2009) Influence of exercise and joint topography on depth-related spatial distribution of proteoglycan and collagen content in immature equine articular cartilage. Equine veterinary journal 41 (6):557-563. doi:10.2746/042516409x424162
Brama PA, TeKoppele JM, Bank RA, Barneveld A, van Weeren PR (2002) Biochemical development of subchondral bone from birth until age eleven months and the influence of physical activity. Equine veterinary journal 34 (2):143-149. doi:10.2746/042516402776767150
Goldring MB, Goldring SR (2010) Articular cartilage and subchondral bone in the pathogenesis of osteoarthritis. Annals of the New York Academy of Sciences 1192:230-237. doi:10.1111/j.1749-6632.2009.05240.x
Ding C, Cicuttini F, Jones G (2007) Tibial subchondral bone size and knee cartilage defects: relevance to knee osteoarthritis. Osteoarthritis and cartilage 15 (5):479-486. doi:10.1016/j.joca.2007.01.003
Kraus VB, Feng S, Wang S, White S, Ainslie M, Brett A, Holmes A, Charles HC (2009) Trabecular morphometry by fractal signature analysis is a novel marker of osteoarthritis progression. Arthritis and rheumatism 60 (12):3711-3722. doi:10.1002/art.25012
Kon E, Filardo G, Brittberg M, Busacca M, Condello V, Engebretsen L, Marlovits S, Niemeyer P, Platzer P, Posthumus M, Verdonk P, Verdonk R, Victor J, van der Merwe W, Widuchowski W, Zorzi C, Marcacci M (2018) A multilayer biomaterial for osteochondral regeneration shows superiority vs microfractures for the treatment of osteochondral lesions in a multicentre randomized trial at 2 years. Knee surgery, sports traumatology, arthroscopy : official journal of the ESSKA 26 (9):2704-2715. doi:10.1007/s00167-017-4707-3
Chen D, Shen J, Zhao W, Wang T, Han L, Hamilton JL, Im HJ (2017) Osteoarthritis: toward a comprehensive understanding of pathological mechanism. Bone research 5:16044. doi:10.1038/boneres.2016.44
Lu W, Li G (2018) Identification of Key Genes and Pathways in Rheumatoid Arthritis Gene Expression Profile by Bioinformatics. Acta reumatologica portuguesa 43 (2):109-131
Zhou J, Li C, Yu A, Jie S, Du X, Liu T, Wang W, Luo Y (2019) Bioinformatics analysis of differentially expressed genes involved in human developmental chondrogenesis. Medicine 98 (27):e16240. doi:10.1097/md.0000000000016240
Wang S, Wang H, Liu W, Wei B (2019) Identification of Key Genes and Pathways Associated with Sex Differences in Osteoarthritis Based on Bioinformatics Analysis. BioMed research international 2019:3482751. doi:10.1155/2019/3482751
He P, Zhang Z, Liao W, Xu D, Fu M, Kang Y (2016) Screening of gene signatures for rheumatoid arthritis and osteoarthritis based on bioinformatics analysis. Molecular medicine reports 14 (2):1587-1593. doi:10.3892/mmr.2016.5423
Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, Marshall KA, Phillippy KH, Sherman PM, Holko M (2012) NCBI GEO: archive for functional genomics data sets—update. Nucleic acids research 41 (D1):D991-D995. doi:10.1093/nar/gks1193
Chou CH, Wu CC, Song IW, Chuang HP, Lu LS, Chang JH, Kuo SY, Lee CH, Wu JY, Chen YT, Kraus VB, Lee MT (2013) Genome-wide expression profiles of subchondral bone in osteoarthritis. Arthritis research & therapy 15 (6):R190. doi:10.1186/ar4380
Korpelainen E, Tuimala J, Somervuo P, Huss M, Wong G (2014) RNA-seq data analysis: a practical approach. CRC press. doi:10.1201/b17457
Smyth GK (2005) Limma: linear models for microarray data. In: Bioinformatics and computational biology solutions using R and Bioconductor. Springer, pp 397-420. doi:10.1007 / 0-387-29362-0_23
Reiner-Benaim A (2007) FDR control by the BH procedure for two-sided correlated tests with implications to gene expression data analysis. Biometrical journal Biometrische Zeitschrift 49 (1):107-126. doi:10.1002/bimj.200510313
Huang DW, Sherman BT, Tan Q, Collins JR, Alvord WG, Roayaei J, Stephens R, Baseler MW, Lane HC, Lempicki RA (2007) The DAVID Gene Functional Classification Tool: a novel biological module-centric algorithm to functionally analyze large gene lists. Genome biology 8 (9):R183. doi:10.1186/gb-2007-8-9-r183
Kanehisa M, Goto S (2000) KEGG: kyoto encyclopedia of genes and genomes. Nucleic acids research 28 (1):27-30. doi:10.1093/nar/28.1.27
Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, Simonovic M, Roth A, Santos A, Tsafou KP, Kuhn M, Bork P, Jensen LJ, von Mering C (2015) STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res 43 (Database issue):D447-452. doi:10.1093/nar/gku1003
Bader GD, Hogue CW (2003) An automated method for finding molecular complexes in large protein interaction networks. BMC bioinformatics 4:2. doi:10.1186/1471-2105-4-2
Wagner AH, Coffman AC, Ainscough BJ, Spies NC, Skidmore ZL, Campbell KM, Krysiak K, Pan D, McMichael JF, Eldred JM, Walker JR, Wilson RK, Mardis ER, Griffith M, Griffith OL (2016) DGIdb 2.0: mining clinically relevant drug-gene interactions. Nucleic Acids Res 44 (D1):D1036-1044. doi:10.1093/nar/gkv1165
Liu W, He J, Lin R, Liang J, Luo Q (2016) Differential proteomics of the synovial membrane between bilateral and unilateral knee osteoarthritis in surgeryinduced rabbit models. Molecular medicine reports 14 (3):2243-2249. doi:10.3892/mmr.2016.5511
Trachana V, Mourmoura E, Papathanasiou I, Tsezou A (2019) Understanding the role of chondrocytes in osteoarthritis: utilizing proteomics. Expert review of proteomics 16 (3):201-213. doi:10.1080/14789450.2019.1571918
Smith T, Baltzer W, Löhr C, Stieger-Vanegas S (2012) Primary synovial osteochondromatosis of the stifle in an English Mastiff. Veterinary and Comparative Orthopaedics and Traumatology 25 (02):160-166
Weber A, Chan PMB, Wen C (2019) Do immune cells lead the way in subchondral bone disturbance in osteoarthritis? Progress in biophysics and molecular biology 148:21-31. doi:10.1016/j.pbiomolbio.2017.12.004
Zhen G, Wen C, Jia X, Li Y, Crane JL, Mears SC, Askin FB, Frassica FJ, Chang W, Yao J, Carrino JA, Cosgarea A, Artemov D, Chen Q, Zhao Z, Zhou X, Riley L, Sponseller P, Wan M, Lu WW, Cao X (2013) Inhibition of TGF-beta signaling in mesenchymal stem cells of subchondral bone attenuates osteoarthritis. Nature medicine 19 (6):704-712. doi:10.1038/nm.3143
Chen Y, Huang YC, Yan CH, Chiu KY, Wei Q, Zhao J, Guo XE, Leung F, Lu WW (2017) Abnormal subchondral bone remodeling and its association with articular cartilage degradation in knees of type 2 diabetes patients. Bone research 5:17034. doi:10.1038/boneres.2017.34
Liang B, Li C, Zhao J (2016) Identification of key pathways and genes in colorectal cancer using bioinformatics analysis. Medical Oncology 33 (10):111
Hepler C, Vishvanath L, Gupta RK (2017) Sorting out adipocyte precursors and their role in physiology and disease. Genes & development 31 (2):127-140. doi:10.1101/gad.293704.116
Pittenger MF, Mackay AM, Beck SC, Jaiswal RK, Douglas R, Mosca JD, Moorman MA, Simonetti DW, Craig S, Marshak DR (1999) Multilineage potential of adult human mesenchymal stem cells. Science (New York, NY) 284 (5411):143-147. doi:10.1126/science.284.5411.143
Li HZ, Xu XH, Lin N, Wang DW, Lin YM, Su ZZ, Lu HD (2020) Overexpression of miR-10a-5p facilitates the progression of osteoarthritis. Aging 12. doi:10.18632/aging.102989
Li Z, Zhong L, Du Z, Chen G, Shang J, Yang Q, Wang Q, Song Y, Zhang G (2019) Network Analyses of Differentially Expressed Genes in Osteoarthritis to Identify Hub Genes. BioMed research international 2019:8340573. doi:10.1155/2019/8340573
Yi CQ, Ma CH, Xie ZP, Cao Y, Zhang GQ, Zhou XK, Liu ZQ (2013) Comparative genome-wide gene expression analysis of rheumatoid arthritis and osteoarthritis. Genetics and molecular research : GMR 12 (3):3136-3145. doi:10.4238/2013.March.11.3
Shirwany NA, Zou MH (2014) AMPK: a cellular metabolic and redox sensor. A minireview. Frontiers in bioscience (Landmark edition) 19:447-474. doi:10.2741/4218
Zheng Z, Xiang S, Wang Y, Dong Y, Li Z, Xiang Y, Bian Y, Feng B, Yang B, Weng X (2020) NR4A1 promotes TNFalphainduced chondrocyte death and migration injury via activating the AMPK/Drp1/mitochondrial fission pathway. International journal of molecular medicine 45 (1):151-161. doi:10.3892/ijmm.2019.4398
Park HY, Kwon HM, Lim HJ, Hong BK, Lee JY, Park BE, Jang Y, Cho SY, Kim HS (2001) Potential role of leptin in angiogenesis: leptin induces endothelial cell proliferation and expression of matrix metalloproteinases in vivo and in vitro. Experimental & molecular medicine 33 (2):95-102. doi:10.1038/emm.2001.17
Yang R, Barouch LA (2007) Leptin signaling and obesity: cardiovascular consequences. Circulation research 101 (6):545-559. doi:10.1161/circresaha.107.156596
Foryst-Ludwig A, Kreissl MC, Benz V, Brix S, Smeir E, Ban Z, Januszewicz E, Salatzki J, Grune J, Schwanstecher AK, Blumrich A, Schirbel A, Klopfleisch R, Rothe M, Blume K, Halle M, Wolfarth B, Kershaw EE, Kintscher U (2015) Adipose Tissue Lipolysis Promotes Exercise-induced Cardiac Hypertrophy Involving the Lipokine C16:1n7-Palmitoleate. The Journal of biological chemistry 290 (39):23603-23615. doi:10.1074/jbc.M115.645341
Strom K, Gundersen TE, Hansson O, Lucas S, Fernandez C, Blomhoff R, Holm C (2009) Hormone-sensitive lipase (HSL) is also a retinyl ester hydrolase: evidence from mice lacking HSL. FASEB journal : official publication of the Federation of American Societies for Experimental Biology 23 (7):2307-2316. doi:10.1096/fj.08-120923
Londos C, Sztalryd C, Tansey JT, Kimmel AR (2005) Role of PAT proteins in lipid metabolism. Biochimie 87 (1):45-49. doi:10.1016/j.biochi.2004.12.010
Francin PJ, Abot A, Guillaume C, Moulin D, Bianchi A, Gegout-Pottie P, Jouzeau JY, Mainard D, Presle N (2014) Association between adiponectin and cartilage degradation in human osteoarthritis. Osteoarthritis and cartilage 22 (3):519-526. doi:10.1016/j.joca.2014.01.002
Cuzdan Coskun N, Ay S, Evcik FD, Oztuna D (2017) Adiponectin: is it a biomarker for assessing the disease severity in knee osteoarthritis patients? International journal of rheumatic diseases 20 (12):1942-1949. doi:10.1111/1756-185x.12790
Gu HY, Yang M, Guo J, Zhang C, Lin LL, Liu Y, Wei RX (2019) Identification of the Biomarkers and Pathological Process of Osteoarthritis: Weighted Gene Co-expression Network Analysis. Frontiers in physiology 10:275. doi:10.3389/fphys.2019.00275

Download PDF

Version 1

posted

You are reading this latest preprint version

Identification of Key Genes and Pathways Associated with Differences of Subchondral Bone in Osteoarthritis

Status:

Version 1

Abstract

Figures

Introduction

Materials And Methods

Results

Discussion

Conclusions

Abbreviations

Declarations

References

Status:

Version 1