Exploring the Key Genes and Pathways of Osteoarthritis in Knee Cartilage in a Rat Model Using Gene Expression Profiling.

PURPOSE: To compare differentially expressed genes (DEGs) mediating osteoarthritis (OA) in knee cartilage and in normal knee cartilage in a rat model of OA and to identify their impact on molecular pathways associated with OA.
MATERIALS AND METHODS: A gene expression profile was downloaded from the Gene Expression Omnibus database. Analysis of DEGs was carried out using GEO2R. Enrichment analyses were performed on the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes pathway using the Search Tool for the Retrieval of Interacting Genes database (http://www.string-db.org/). Subsequently, the regulatory interaction network of OA-associated genes was visualized using Cytoscape software (version 3.4.0; www.cytoscape.org).
RESULTS: In the gene expression profile GSE103416, a total of 99 DEGs were identified. Among them, 76 DEGs (76.77%) were overexpressed, and the remaining 23 DEGs (23.23%) were underexpressed. GO and pathway enrichment analyses of target genes were performed. Using gene-gene interaction network analysis, relevant core genes, including MET, UBB, GNAI3, and GNA13, were shown to hold a potential relationship with the development of OA in cartilage. Using quantitative real-time PCR, the Gna13/cGMP-PKG signaling pathway was identified as a potential research target for therapy and for further understanding the development of OA.
CONCLUSION: The results of the present study provide a comprehensive understanding of the roles of DEGs in knee cartilage in relation to the development of OA. © Copyright: Yonsei University College of Medicine 2018.

INTRODUCTION

Knee osteoarthritis (OA) is a chronic degenerative disease of the joints that is common in older adults, and its incidence increases with aging.1 Reportedly, the total prevalence rate of OA is 15%, with rates of 10–17% in patients over 40 years, 50% in patients over 60 years, and as much as 80% in patients 75 years of age.2 OA can lead to dysfunction and knee pain. Although it is not fatal, OA seriously affects patient quality of life. The pathogenesis of OA is complex, involving cartilaginous aging, bone remodeling, synovial lesions, microvascular lesions, inflammation, and other biological processes.345 Its specific pathological mechanisms, however, have not yet been elucidated.

The pathogenesis of OA mainly involves pathological changes in articular cartilage and chondrocytes. In recent years, molecular biology has played an important role in explaining its disease pathophysiology. Gene regulation has been implicated in driving an imbalance between the expression of catabolic and anabolic factors, leading eventually to osteoarthritic cartilage degeneration. The various stages of disease progression have been described by a complex pattern of transcriptional regulations.6

In the present study, we downloaded the gene expression profile GSE103416 from the Gene Expression Omnibus database (GEO) and conducted bioinformatics analysis to investigate differentially expressed genes (DEGs) between OA knee cartilage and normal knee cartilage in a rat model of OA. In doing so, we hoped that the key target genes and pathways involved in the progression of OA could be identified and that existing molecular mechanisms could be revealed.

MATERIALS AND METHODS

Gene expression microarray data

Knee cartilage samples were collected from individual male Wistar rats with monosodium iodoacetate (MIA)-induced OA (0, 2, 14, and 28 days after treatment). MIA-treated (3 mg) animals were compared to intact controls. Four replicate samples were analyzed per experimental group. The gene expression profile GSE103416 was downloaded from the GEO (www.ncbi.nlm.nih.gov/geo/). GSE103416 was based on the Affymetrix Rat Gene 2.0 ST Array GPL17117 platform. The GSE103416 dataset contained 16 samples, including 12 OA knee cartilage samples and four normal control knee cartilage samples. A box plot of the sample data is provided in Fig. 1.

Fig. 1

Box plot for the sample data after normalization.

DEGs in OA knee cartilage and control cartilage

The raw data files used for the analysis included TXT files (Affymetrix platform). The analysis was carried out using GEO2R, which can perform comparisons on original submitter-supplied processed data tables using the GEO query and limma R packages from the Bioconductor project. A p-value <0.05, adjusted (adj.) p-value <0.05, and log fold change (FC) >1.0 and log FC <−1.0 were used as cut-off criteria. DEGs with statistical significance between OA knee cartilage and control cartilage were selected and identified.

GO and the KEGG analysis of DEGs

The target gene list was subjected to analysis using Cytoscape software (version 3.4.0) (www.cytoscape.org) and ClueGO (version 2.33, Laboratory of Integrative Cancer Immunology, Bethesda, MD, USA) to identify overrepresented Gene Ontology (GO) categories and pathway categories. GO analysis was used to predict the potential functions of the DEGs in biological processes (BP), molecular functions (MF), and cellular components (CC). Finally, the overrepresented pathway categories with p-value <0.05 were considered statistically significant using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis: KEGG (http://www.genome.jp/) is a knowledge base for systematic analysis of gene functions, linking genomic information with higher-level systemic functions.

Gene interaction network construction

A large number of the DEGs we obtained were potential OA-related-genes, and it was suggested that these DEGs in knee cartilage may participate in the progression of OA. First, the DEGs list was subjected to analysis by the Search Tool for the Retrieval of Interacting Genes (STRING) database (http://www.string-db.org/), and an interaction network chart with a combined score >0.4 was saved and exported. Subsequently, the interaction regulatory network of OA-associated genes was visualized using Cytoscape software (version 3.4.0). Then, the plugin Molecular Complex Detection (MCODE) was applied to screen the modules of the gene interaction network in Cytoscape.

qRT-PCR verification

Six genes were amplified by quantitative real-time polymerase chain reaction (qRT-PCR) with the Bio-Rad CFX96 Real-Time PCR System (Bio-Rad, Invitrogen Corp., Carlsbad, CA, USA). Total RNA was isolated from the MIA-induced OA knee cartilage samples and normal knee cartilage samples using TRIzol reagent (Invitrogen Corp, Carlsbad, CA, USA), and was polyadenylated and reverse-transcribed with a poly(T) adapter into cDNA, following the manufacturer's directions. Real-time PCR was performed using SYBR green dye in a thermal cycler with the following parameters: an initial denaturation step at 95℃ for 30 min; 40 cycles at 95℃ for 5 seconds; and 60℃ for 30 seconds. The complete experimental process was performed for each sample in triplicate. All primers were synthesized by Shanghai Shenggong Inc. (Shanghai, China) and mRNA-specific primers are listed in Supplementary Table 1 (only online). All data were analyzed using the 2−ΔΔCT method to calculate the difference between the threshold cycle (CT) values of the target genes in each sample. All data were analyzed using SPSS statistical software (version 11.5 for Windows, SPSS Inc., Chicago, IL, USA). Statistical analysis was performed using two-tailed Student's t-test, and differences were considered statistically significant at p<0.05.

RESULTS

Identification of DEGs

The gene expression profile GSE103416 was downloaded from GEO, and the GEO2R method was used to compare the DEGs in the OA knee cartilage samples with those of the intact controls. In the OA knee cartilage samples, 99 DEGs were identified, including 76 upregulated DEGs and 23 downregulated DEGs, compared with normal knee cartilage samples. The volcano plot and partial list DEGs are shown in Fig. 2 and Table 1, respectively.

Fig. 2

The volcano plot of DEGs. The abscissa indicates log2 FC, and the ordinates are −log10 (p-value). Each point represents a gene. The green dots represent the downregulated DEGs, red dots represent the upregulated DEGs, and black dots represent non-DEGs. DEG, differentially expressed gene; FC, fold change.

Table 1

The Top 10 Regulated Differentially Expressed Genes between Osteoarthritis Knee Cartilage and the Normal Knee Cartilage

adj. p value	p value	log FC	Gene symbol
Upregulated
1.51E-05	3.49E-08	3.10903306	CILP
6.23E-04	9.47E-06	2.71209289	MYOC
3.77E-06	2.26E-09	2.66375242	VSIG4
4.53E-05	1.95E-07	2.52905644	WIF1
2.26E-03	6.70E-05	2.48010086	SCD
6.51E-06	6.91E-09	2.42215222	ELOVL7
2.13E-02	1.84E-03	2.39852831	COL2A1
3.21E-03	1.13E-04	2.29479274	CHAD
1.87E-03	5.03E-05	2.20344792	AOX3
6.42E-03	3.09E-04	2.19130138	PLIN1
Downregulated
1.52E-03	3.79E-05	−1.64626895	CX3CR1
2.46E-02	2.27E-03	−1.45059698	Clec4a3
9.75E-03	5.72E-04	−1.42840862	P4HA3
3.26E-02	3.41E-03	−1.36439436	CXCL6
2.66E-05	8.17E-08	−1.3377121	PTPRN
1.40E-02	9.99E-04	−1.3150188	PHYHIPL
3.36E-02	3.56E-03	−1.30002041	CD8B
4.31E-05	1.78E-07	−1.28838008	TFPI2
6.76E-05	3.56E-07	−1.2799337	GPNMB
9.89E-03	5.87E-04	−1.22216962	KCNJ15

FC, fold change.

GO term enrichment analysis of DEGs

Functional annotation of the 99 DEGs was clarified using the Cytoscape software online tool. GO analysis (Fig. 3) indicated that these DEGs were significantly enriched in protein metabolism, signal transduction, cell communication, transport, energy pathways, metabolism, and other BP. For MF, the DEGs were enriched in activity for glucosidase, oxidoreductase, B cell receptor, and translation regulator activities, as well as activities of other receptors. In addition, GO CC analysis showed that the DEGs were significantly enriched in the Golgi apparatus, integral to the plasma membrane, lysosomal membranes, endoplasmic reticulum, and others.

Fig. 3

Gene Ontology-enrichment analysis of biological processes (A), molecular functions (B), and cellular components (C). OA, osteoarthritis.

KEGG pathway analysis of DEGs

The result of KEGG pathway analysis revealed that target genes were enriched in steroid biosynthesis, regulation of lipolysis in adipocytes, renin secretion, apoptosis, ether lipid metabolism, the TNF signaling pathway, endocrine and other factor-regulated calcium reabsorption, protein processing in the endoplasmic reticulum, the cGMP-PKG signaling pathway, thyroid hormone synthesis, cholinergic synapses, the adipocytokine signaling pathway, the NF-kappa B signaling pathway, the estrogen signaling pathway, salivary secretion, glycerophospholipid metabolism, the thyroid hormone signaling pathway, the cAMP signaling pathway, the phospholipase D signaling pathway, regulation of lipolysis in adipocytes, the PI3K-Akt signaling pathway, and the MAPK signaling pathway. Most of the pathways were not reported in any previous study, and only the NF-kappa B signaling pathway, the PI3K-Akt signaling pathway, and the MAPK signaling pathway would be considered well-known pathways in the development or progression of OA. These key pathways are shown in Fig. 4. These core pathways and their associated genes are summarized in Table 2. The first-ranking steroid biosynthesis signaling pathway accounted for 35% of the associated genes, including CYP51A1, DHCR24, EBP, LSS, MSMO1, SOAT1, and SQLE.

Fig. 4

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of differentially expressed genes (DEGs). The different node colors mean different pathways, and the closer the colors are, the closer the function clustering of pathways are.

Table 2

Core Pathways and Their Associated Genes Found

GO ID	GO Term	Term p value	% associated genes	Associated genes found
GO:0000100	Steroid biosynthesis	4.50E-05	35.00	[CYP51A1, DHCR24, EBP, LSS, MSMO1, SOAT1, SQLE]
GO:0004923	Regulation of lipolysis inadipocytes	9.20E-05	20.37	[ADCY5, ADCY6, ADORA1, ADRB3, AKT1, AKT2, FABP4, PLIN1, PTGS1, PTGS2, TSHR]
GO:0004924	Renin secretion	2.60E-05	20.00	[ADCY5, ADCY6, ADCYAP1R1, ADORA1, ADRB3, CACNA1D, GUCY1A2, GUCY1A3, ITPR2, KCNJ2, PLCB4, PTGER2, PTGER4]
GO:0004215	Apoptosis	6.70E-03	18.18	[BAK1, BCL2, BCL2L11, CASP3, PMAIP1, TNFRSF1A]
GO:0000565	Ether lipid metabolism	2.90E-02	13.33	[CHPT1, ENPP2, LPCAT2, PAFAH1B3, PLD1, PLPP1]
GO:0004668	TNF signaling pathway	1.50E-03	12.96	[AKT1, AKT2, CASP3, CCL20, CREB3, CX3CL1, ITCH, MAP2K4, MAP3K5, NFKB1, PTGS2, SELE, TNFRSF1A, TNFRSF1B]
GO:0004961	Endocrine and other factor-regulated calcium reabsorption	3.50E-02	12.77	[ADCY6, ATP1B1, DNM3, ESR1, KL, PLCB4]
GO:0004141	Protein processing in endoplasmic reticulum	1.50E-04	12.65	[BAK1, BCL2, CKAP4, DNAJC3, EIF2S1, ERO1B, FBXO6, HSP90AB1, HSP90B1, HSPBP1, LMAN1, MAP3K5, MOGS, PDIA4, PRKCSH, PRKN, RPN2, RRBP1, TXNDC5, UBE4B, YOD1]
GO:0004022	cGMP-PKG signaling pathway	4.60E-04	12.27	[ADCY5, ADCY6, ADORA1, ADRA1A, ADRA1D, ADRB3, AKT1, AKT2, ATP1B1, ATP2A3, CACNA1D, CREB3, GNA13, GUCY1A2, GUCY1A3, ITPR2, PLCB4, PLN, PPP1R12A, SLC25A5]
GO:0004918	Thyroid hormone synthesis	1.50E-02	12.16	[ADCY5, ADCY6, ATP1B1, CREB3, HSP90B1, ITPR2, PDIA4, PLCB4, TSHR]
GO:0004725	Cholinergic synapse	8.00E-03	11.61	[ADCY5, ADCY6, AKT1, AKT2, BCL2, CACNA1D, CHRM5, CHRNB4, CREB3, ITPR2, JAK2, KCNJ2, PLCB4]
GO:0004920	Adipocytokine signaling pathway	2.70E-02	11.59	[AKT1, AKT2, JAK2, NFKB1, PRKAG2, SLC2A4, TNFRSF1A, TNFRSF1B]
GO:0004064	NF-kappa B signaling pathway	1.60E-02	11.58	[BCL2, BCL2A1, CD14, LY96, MYD88, NFKB1, PLAU, PTGS2, TLR4, TNFRSF13C, TNFRSF1A]
GO:0004915	Estrogen signaling pathway	1.80E-02	11.22	[ADCY5, ADCY6, AKT1, AKT2, CREB3, ESR1, HBEGF, HSP90AB1, HSP90B1, ITPR2, PLCB4]
GO:0004970	Salivary secretion	2.60E-02	11.11	[ADCY5, ADCY6, ADRA1A, ADRA1D, ADRB3, ATP1B1, GUCY1A2, GUCY1A3, ITPR2, PLCB4]
GO:0000564	Glycerophospholipid metabolism	3.20E-02	10.53	[CHPT1, LPCAT2, LYPLA1, LYPLA2, MBOAT7, PCYT1A, PGS1, PLD1, PLPP1, PTDSS1]
GO:0004919	Thyroid hormone signaling pathway	2.00E-02	10.34	[AKT1, AKT2, ATP1B1, ESR1, HDAC2, MED12L, NCOA2, NCOR1, PLCB4, PLN, SLC16A10, THRB]
GO:0004024	cAMP signaling pathway	5.10E-03	10.10	[ABCC4, ADCY5, ADCY6, ADCYAP1R1, ADORA1, AFDN, AKT1, AKT2, ATP1B1, CACNA1D, CREB3, GRIA2, NFKB1, PLD1, PLN, PPP1R12A, PTGER2, RAPGEF3, RYR2, TSHR]
GO:0004072	Phospholipase D signaling pathway	4.E-02	9.59	[ADCY5, ADCY6, AKT1, AKT2, DNM3, GNA13, GRM8, KIT, PLCB4, PLD1, PLPP1, RALB, RALGDS, RAPGEF3]
GO:0004923	Regulation of lipolysis in adipocytes	9.20E-05	20.37	[ADCY5, ADCY6, ADORA1, ADRB3, AKT1, AKT2, FABP4, PLIN1, PTGS1, PTGS2, TSHR]

Interaction network of DEGs and core genes in the interaction network

Based on information in the STRING database, the gene interaction network contained 757 nodes and 1947 edges. The nodes indicated the DEGs, and the edges indicated the interactions among the DEGs. NetworkAnalyzer in Cytoscape software was used to analyze these genes, and core genes were ranked according to the predicted scores. The top 20 high-degree hub nodes included AKT1, HSP90AB1, PHLPP1, CASP3, BCL2, SMARCA4, HDAC2, JAK2, GNAI3, PIK3CB, PTGS2, ESR1, TLR4, AKT2, CDK4, MET, UBB, ITGAM, GNA13, and NFKB1. Among these, AKT1, CASP3, BCL2, JAK2, PTGS2, ESR1, TLR4, and AKT2 are well known in the development or progression of OA. AKT1 showed the highest node degree of 86. The core genes and their corresponding degrees are shown in Table 3. The distribution of core genes in the interaction network is shown in Fig. 5. The correlation between data points and corresponding points on the line was approximately 0.893. The R-squared value was 0.887, suggesting relatively high confidence that the underlying model was indeed linear. Then, we used MCODE to screen the modules of the gene interaction network, and the top three modules are shown in Fig. 6.

Table 3

The Core Genes and Their Corresponding Degree

Gene	Degree	Gene	Degree	Gene	Degree	Gene	Degree
Akt1	86	Smarca4	32	Ptgs2	36	Met	30
Hsp90ab1	66	Hdac2	31	Esr1	35	Ubb	30
Phlpp1	59	Jak2	29	Tlr4	34	Itgam	28
Casp3	45	Gnai3	28	Akt2	34	Gna13	27
Bcl2	44	Pik3cb	37	Cdk4	31	Nfkb1	27

Fig. 5

The distribution of core genes in the interaction network. The black node means the core gene. The red line means the fitted line, and the blue line means the power law. The correlation between the data points and corresponding points on the line is approximately 0.893. The R-squared value is 0.887, giving a relatively high confidence that the underlying model is indeed linear.

Fig. 6

The top three modules from the gene-gene interaction network. The squares represent the differentially expressed genes (DEGs) in modules, and the lines show the interaction between the DEGs.

Genes expression validation by qRT-PCR

To validate the bioinformatic analysis results, qRT-PCR was used to quantify parts of the explored genes, including two upregulated genes (MET, UBB) and two downregulated genes (GNAI3, GNA13). As shown in Fig. 7, the mRNA expression patterns of Met, Ubb, Gnai3, and Gna13 detected by qRT-PCR significantly accorded with the corresponding gene alteration of microarray data (p<0.05). In addition, mRNA expression levels of PKGII and Sox9 in OA knee cartilage were detected by qRT-PCR: PKGII was downregulated and Sox9 was upregulated in the OA knee cartilage groups, compared with the normal knee cartilage groups (p<0.05).

Fig. 7

Validation of the differential expression of mRNA of corresponding genes identified in the OA knee cartilage groups, compared with NC groups, by qRT-PCR. Data indicate relative expression following normalization. Values are means±standard error (*p<0.05). Figures (A–C) show that three mRNAs of differentially expressed genes (Gnai3, Gna13, PKGII) were downregulated in OA knee cartilage groups compared with normal knee cartilage groups. Figures (D–F) showed that three mRNAs of genes (Met, Ubb, Sox9) were upregulated in OA knee cartilage groups, compared with normal knee cartilage groups. OA, osteoarthritis; NC, normal control.

DISCUSSION

OA is a progressive disease of the joints characterized by the breakdown of the cartilage matrix, chondrocyte hypertrophy, and ectopic growth of bony structures in the joints.78 Chondrocyte apoptosis is largely responsible for cartilage degeneration in OA. In the present study, the gene expression profile of GSE103416 was downloaded and a bioinformatics analysis was performed. The results revealed 99 DEGs in OA knee cartilage samples, compared with normal controls. Furthermore, GO, KEGG pathway, and gene-gene interaction network analysis were performed to identify biomarkers or the major genes related to cytogenetic pathways in cartilage degeneration associated with OA.

To elucidate the underlying differences in molecular mechanisms between OA knee cartilage and normal cartilage, we characterized the possible GO functional terms and signaling pathways of DEGs. Considering the results of GO function analysis, we linked the DEGs with the following processes: protein metabolism, signal transduction, and apoptosis, all of which would likely be important for the development of OA. Similar to what previous articles reported, our KEGG pathway analysis showed that the PI3K-Akt signaling pathway, the NF-kappa B signaling pathway, the MAPK signaling pathway, and others were among the most relevant pathways for OA. Xu, et al.9 found that high expression of Sam68 promoted NF-kappa B signaling activation, catabolic gene expression, and cellular apoptosis in TNF-α-treated chondrocytes, providing insight into the pathophysiology of OA, as well as a potential target for its treatment. He and Cheng10 showed that miR-20 had an important role in the pathogenesis of OA, inhibiting chondrocyte proliferation and autophagy by targeting ATG10 via the PI3K/AKT/mTOR signaling pathway. Li, et al.11 found that low-intensity pulsed ultrasound treatment promoted cartilage repair in OA via downregulation of the MAPK signaling pathway. All these signaling pathways have been reported to play important roles in molecular mechanism of degenerative processes in OA knee cartilage.

Notably, there was abundant evidence of DEGs in OA knee cartilage that have been proven to play important roles during the pathological process of OA. The STRING database revealed the top 20 high-degree hub nodes of DEGs, including AKT1, HSP90AB1, PHLPP1, CASP3, BCL2, SMARCA4, HDAC2, JAK2, GNAI3, PIK3CB, PTGS2, ESR1, TLR4, AKT2, CDK4, MET, UBB, ITGAM, GNA13, and NFKB1. Furthermore, we analyzed the gene interaction network and the top three modules using MCODE, and found that PTGS2, ESR1, TLR4, and AKT2 were the core interaction genes for OA. This finding was in accordance with our STRING database results. Chondrocyte apoptosis is largely responsible for cartilage degeneration in OA. Xue, et al.12 found that lactoferrin triggered CREB1 phosphorylation in IL-1β-induced human articular chondrocyte apoptosis through AKT1 signaling. Bradley, et al.13 showed that Phlpp1 deficiency increased Akt2 activity, diminishing FoxO1 levels and inducing Fgf18 expression to stimulate Erk1/2 activity and chondrocyte proliferation. They also found that Phlpp1 deficiency protected against OA progression, while CpG demethylation and inflammatory responses promoted Phlpp1 expression.14 By histological analysis, TUNEL assay and immunohistochemical analyses of Casp3, Bcl-2, and Bax expression, Zhou, et al.15 found that berberine ameliorated cartilage degeneration and exhibited an anti-apoptotic effect on articular cartilage in a rat OA model. In Culley, et al.'s study,16 inhibition of class I HDACs (HDAC-1, HDAC-2, HDAC-3) by MS-275 or by specific depletion of HDACs repressed cytokine-induced metalloproteinase expression in cartilage cells and bovine nasal cartilage explants, resulting in inhibition of cartilage resorption. This result indicated that specific inhibition of class I HDACs was a possible therapeutic strategy in OA. In addition, Zhang, et al.17 found that leptin induced apoptosis of chondrocytes in an in vitro model of OA via the JAK2/STAT3 signaling pathway. Valdes, et al.18 performed a case-control study of men and women with clinical knee OA genetic associations in 12 candidate genes previously reported to be associated with OA, and found that ESR1 and other genes might be helpful in identifying women who were at high risk of developing knee OA. The expression of TLR4 in cartilage was increased throughout the development of OA, and many TLR4 agonists that have been identified in the joints of patients with OA induced inflammatory responses in ex vivo tissue samples from these patients.19 Schneider, et al.20 found that the promoter variant of the COX-2/PTGS2 gene was associated with a lower risk for end-stage hip and knee OA. CDK4 also correlated with the disease grades of knee OA in Zhu, et al.'s studies.21 Taken together, AKT1, PHLPP1, CASP3, BCL2, HDAC2, JAK2, PTGS2, ESR1, TLR4, and AKT2 genes have been reported to take part in the development of OA.

Similar bioinformatics analysis in mouse models or human cartilage has contributed to further revealing genes important in OA. Using GSE48556 data, Feng and Lian22 found that DEGs, such as PDGFRB, IFNG, EGR1, FASLG, and H3F3B, may be potential targets for OA diagnosis and treatment. Using four gene expression profiling datasets and verification of qRT-PCR, Tian23 found that TWIST1, POU2F1, SMARCA4, and CREBBP may play important roles in OA pathology, and SMARCA4 was also reported in our study. Mi, et al.24 reported that their top 10 hub genes were UBC, RPL23A, mammalian target of rapamycin, HSP90AB1, RPS28, RPL37A, RPS24, RPS4X, RPS18, and UBB. Of these genes, HSP90AB1 and UBB were also reported in our study. In addition, another four bioinformatics analyses identified AGT, CXCL12, KDM2B, VEGFA, JUN, JUNB, PISD, RARRES3, EIF4G1, EPHA3, KIF2C, KIF11, KIF20A, PTTG1, MAP2K6, PPP3CC, and CSNK1E as potentially essential genes associated with the pathogenesis of knee OA.25262728

Bioinformatics analysis alone is insufficient, and corresponding in vitro or in vivo experiments should be added.29 To further validate our results, qRT-PCR was performed for two upregulated genes (MET, UBB) and two downregulated (GNAI3, GNA13), which proved expression changes of these four genes were in accordance with the microarray data. Interestingly, according to the KEGG pathway results, GNA13 was associated with the cGMP-PKG signaling pathway, which has not been reported in research on OA development or progression. Previous reports showed that the cGMP-PKG signaling pathway is related to chondrocyte function, shedding light on a novel role of PKGII as a molecular switch, coupling the cessation of proliferation and the start of hypertrophic differentiation of chondrocytes through the attenuation of Sox9 function.3031 Hence, we detected mRNA expression levels of cGMP-dependent protein kinases type II (PKGII) and Sox9 in OA knee cartilage, and verified decreased PKGII expression and increased Sox9 expression in the OA knee cartilage groups, indicating that the cGMP-PKG signaling pathway may contribute to the development of OA in cartilage. Further validation, however, is required.

In summary, 99 DEGs were identified, including 76 upregulated DEGs and 23 downregulated DEGs, between OA knee cartilage and normal cartilage in a rat model. GO and KEGG pathway analysis provided a series of related key genes and pathways to contribute to the understanding of the molecular mechanisms of OA knee cartilage degeneration or apoptosis, yielding clues suggesting that the genes MET, UBB, GNAI3 and GNA13, as well as the cGMP-PKG signaling pathway, are highly correlated with the development of OA in cartilage.