Identification of Pivotal Genes and Pathways in Osteoarthritic Degenerative Meniscal Lesions via Bioinformatics Analysis of the GSE52042 Dataset.

BACKGROUND To better understand the process of osteoarthritic degenerative meniscal lesions (DMLs) formation, this study analyzed the dataset GSE52042 using bioinformatics methods to identify the pivotal genes and pathways related to osteoarthritic DMLs. MATERIAL AND METHODS The GSE52042 dataset, comprising diseased meniscus samples and healthier meniscus samples, was downloaded and the differentially-expressed genes (DEGs) were extracted. The reactome pathways assessment and functional analysis were performed using the "ClusterProfiler" package and "ReactomePA" package of Bioconductor. The protein-protein interaction network was constructed, followed by the extraction of hub genes and modules. RESULTS A set of 154 common DEGs, including 64 upregulated DEGs and 90 downregulated DEGs, were obtained. GO analysis suggested that the DEGs primarily participated in positive regulation of the mitotic cell cycle and extracellular matrix organization. Reactome pathway analysis showed that the DEGs were predominantly enriched in TP53, which regulates transcription of genes involved in G2 cell cycle arrest and extracellular matrix organization. The top 10 hub genes were TYMS, AURKA, CENPN, NUSAP1, CENPM, TPX2, CDK1, UBE2C, BIRC5, and CCNB1. The genes in the 2 modules were primarily associated with M Phase and keratan sulfate degradation. CONCLUSIONS A series of pivotal genes and reactome pathways were identified elucidate the molecular mechanisms involved in the formation of osteoarthritic DMLs and to discover potential therapeutic targets.

Background

The meniscus is an important structure of the knee, playing a pivotal role in shock absorption, joint lubrication and nutrition, proprioception, load transmission, joint congruity, and stability [1]. A meniscus lesion, predisposing to knee osteoarthritis (OA), is a frequent knee joint injury with high incidence rate, found in 60% of patients with asymptomatic osteoarthritis [2-4]. There is a high interdependency between degenerative meniscal lesions (DMLs) and knee OA [5,6]. DMLs can be a feature of knee OA and are conducive to joint space narrowing [7-9]. The best treatment of DMLs is still debated. Although arthroscopic surgery is still widely performed in the management of DMLs [10], guidelines and research recommend against surgery due to the non-superiority of arthroscopic surgery over non-operative treatment [11-13]. Biological therapies such as platelet-rich plasma intrameniscal injection may be promising treatments for DMLs [14]. Therefore, better understanding of biological properties and molecular mechanisms in progression of DMLs may promote development of more effective therapeutic strategies.

Recently, microarray technology, which is an effective and high-throughput molecular technique for simultaneously analyzing expression of multiple genes, has been pervasively used in a range of research areas to discover the differentially-expressed genes within molecular mechanisms, distinct pathways, and protein–protein interactions (PPI), and to reveal connections between disease and genes [15-20]. Sun et al. detected many differentially-expressed genes (DEGs) between osteoarthritic meniscal cells and normal meniscal cells [21]. Brophy et al. found the transcripts and biological processes in OA menisci were different from non-OA menisci [22]. Rai et al. revealed a correlation of body mass index to gene expression in ruptured menisci and a biologic association between obesity and knee osteoarthritis [23]. Microarray analysis combined with bioinformatics provides a novel and effective method to investigate the comprehensive molecular mechanisms in meniscus injuries [24], and may be useful to improve understanding of the connections between DEGs and DMLs.

To explore the molecular mechanisms of DMLs in OA patients, the present study aimed to identify the key DEGs and functional pathways in OA DMLs by use of comprehensive bioinformatics methods. As OA menisci have unique biological properties, we reanalyzed the gene expression profile of GSE52042, which identified the DEGs between diseased or severe degenerative tear OA menisci and healthier OA menisci. Subsequently, extensive bioinformatics analysis was carried out for gene ontology (GO) analysis, reactome pathways identification, PPI network screening, and modules extraction. The results of this work may provide further insights into the pathogenesis of OA DMLs at the molecular level and help explore potential biotherapeutic targets for it. To the best of our knowledge, there is no previously published bioinformatics analysis of this database.

Material and Methods

Microarray data

The GSE52042 dataset in Gene Expression Omnibus () was downloaded for subsequent analysis. GSE52042 was based on GPL17882 platform (Human MI ReadyArray™ 49K Human Genomic Array; Microarrays, Inc., Huntsville, AL, USA), including 4 diseased meniscus samples with thoroughly degraded surfaces and 4 healthier meniscus samples with integrated surfaces, which were obtained from late-stage of knee OA.

Identification of DEGs in diseased and healthier samples

The DEGs between diseased and healthier samples were screened by GEO2R (). P value <0.05 and log fold change (FC) <−1.0 or log FC >1.0 were defined as the cutoff criteria. A volcano plot and a heatmap were constructed using R software (version 3.6.0).

Functional and reactome pathway analysis of DEGs

GO enrichment analysis was carried out by the “ClusterProfiler” package of Bioconductor [25] to assess the potential functions of DEGs in molecular function (MF), cellular component (CC), and biological process (BP). Overrepresented enriched reactome pathways were searched for using the “ReactomePA” package of Bioconductor [26] to identify pathways that were strongly related to DEGs. The P value <0.05 and false discovery rate (FDR) <0.05 were set as the cutoff for criteria.

Construction of PPI network and modules analysis

All DEGs were uploaded to STRING (version 11.0; ) to evaluate the PPI interactive relationships among DEGs. Confidence score >0.4 was defined as the cutoff criterion. Subsequently, the PPI network was visualized using Cytoscape software (version 3.6.1; www.cytoscape.org). The cytoHubba plugin (version 0.1) [27] was used to extract the hub genes. Then, the modules were extracted using the Molecular Complex Detection plugin (version 1.5.1) according to the criteria: Degree Cutoff=10, Max. Depth=100, K-Core=2, Mode Score Cutoff=0.2, and MCODE score ≥4 [28].

Results

Identification of DEGs

In comparing diseased samples to healthier samples, a total of 154 DEGs including 64 upregulated and 90 downregulated genes were identified from the GSE52042 dataset. The top 10 upregulated and downregulated DEGs are shown in Table 1. The heatmap and volcano plot for total DEGs are illustrated in Figures 1 and 2, respectively.

Table 1

The top 10 upregulated and downregulated DEGs.

	Gene symbol	logFC	P value
Downregulated	FMOD	−2.396	0.000533
	MXRA5	−2.255	0.021918
	VCAM1	−1.923	0.000582
	LRP1	−1.922	0.011987
	CFH	−1.879	0.002062
	TIMP2	−1.878	5.67E-05
	OAS1	−1.844	0.000157
	CSPG4	−1.835	0.000205
	LTBP3	−1.804	0.002253
	THBS2	−1.726	0.007495
Upregulated	PSG9	2.081	0.003259
	PSG4	1.927	0.003893
	PSG11	1.836	0.004486
	FSTL1	1.776	0.001282
	PSG8	1.774	0.000716
	EGR1	1.739	9.56E-05
	FOS	1.732	0.000996
	NEK7	1.728	0.015207
	ACTR2	1.683	0.004529
	KRT18P21	1.608	7.42E-06

Figure 1

Heat map plot of total DEGs.

Figure 2

Volcano plot of total DEGs.

GO terms enrichment analysis

In terms of BP, the upregulated DEGs were primarily participated in female pregnancy, positive regulation of cell cycle, multi-multicellular organism process, positive regulation of mitotic cell cycle, and chromosome segregation, and the downregulated DEGs were chiefly enriched in extracellular matrix (ECM) organization, extracellular structure organization, glycosaminoglycan catabolic process, aminoglycan catabolic process, and keratan sulfate catabolic process.

In terms of CC, the upregulated DEGs were mostly enriched in spindle, midbody, spindle pole, mitotic spindle, and condensed chromosome kinetochore, and the downregulated DEGs were predominantly enriched in collagen-containing extracellular matrix, extracellular matrix, lysosomal lumen, vacuolar lumen, and Golgi lumen.

In terms of MF, the upregulated DEGs were predominantly enriched in ubiquitin-like protein ligase binding. The downregulated DEGs were chiefly involved in extracellular matrix structural constituent, extracellular matrix structural constituent conferring compression resistance, growth factor binding, collagen binding, and glycosaminoglycan binding.

The top 5 BP, CC, and MF enrichment analyses of the DEGs are displayed in Table 2. The bubble charts of the DEGs are shown in Figures 3 and 4, which exhibit the overview of GO terms (the top 5 BP, CC, and MF enrichment analyses). The chord diagrams of the DEGs are illustrated in Figures 5 and 6, which display the relationship between DEGs and GO terms (the top 5 BP, CC, and MF enrichment analyses).

Table 2

The top 5 BP, CC, and MF enrichment analyses of the upregulated and downregulated DEGs.

Ontology	GO ID	Description	Count	Gene ratio	P value	Adjust p value	FDR
Downregulated
BP	GO: 0030198	Extracellular matrix organization	18	0.22	9.33E-15	2.08E-11	1.57E-11
BP	GO: 0043062	Extracellular structure organization	18	0.22	1.18E-13	1.31E-10	9.87E-11
BP	GO: 0006027	Glycosaminoglycan catabolic process	8	0.10	3.29E-10	2.39E-07	1.80E-07
BP	GO: 0006026	Aminoglycan catabolic process	8	0.10	4.29E-10	2.39E-07	1.80E-07
BP	GO: 0042340	Keratan sulfate catabolic process	5	0.06	1.37E-09	4.85E-07	3.66E-07
CC	GO: 0062023	Collagen-containing extracellular matrix	24	0.30	1.59E-21	2.83E-19	2.17E-19
CC	GO: 0031012	Extracellular matrix	24	0.30	6.65E-20	5.92E-18	4.55E-18
CC	GO: 0043202	Lysosomal lumen	10	0.13	6.91E-12	4.10E-10	3.15E-10
CC	GO: 0005775	Vacuolar lumen	11	0.14	1.45E-10	6.43E-09	4.95E-09
CC	GO: 0005796	Golgi lumen	8	0.10	1.13E-08	4.04E-07	3.10E-07
MF	GO: 0005201	Extracellular matrix structural constituent	15	0.19	6.88E-16	1.64E-13	1.36E-13
MF	GO: 0030021	Extracellular matrix structural constituent conferring compression resistance	5	0.06	3.94E-08	4.68E-06	3.89E-06
MF	GO: 0019838	Growth factor binding	7	0.09	1.63E-06	0.000129	0.000107
MF	GO: 0005518	Collagen binding	4	0.05	0.000195	0.011595	0.009641
MF	GO: 0005539	Glycosaminoglycan binding	6	0.08	0.000353	0.016785	0.013956
Upregulated
BP	GO: 0007565	Female pregnancy	8	0.14	1.43E-07	0.000172	0.000127
BP	GO: 0045787	Positive regulation of cell cycle	10	0.18	2.06E-07	0.000172	0.000127
BP	GO: 0044706	Multi-multicellular organism process	8	0.14	4.63E-07	0.000206	0.000151
BP	GO: 0045931	Positive regulation of mitotic cell cycle	7	0.13	4.93E-07	0.000206	0.000151
BP	GO: 0007059	Chromosome segregation	8	0.14	2.09E-06	0.000699	0.000513
CC	GO: 0005819	Spindle	8	0.14	7.09E-06	0.001305	0.001105
CC	GO: 0030496	Midbody	6	0.10	1.47E-05	0.001354	0.001146
CC	GO: 0000922	Spindle pole	5	0.09	0.000108	0.006634	0.005617
CC	GO: 0072686	Mitotic spindle	4	0.07	0.000234	0.009674	0.008191
CC	GO: 0000777	Condensed chromosome kinetochore	4	0.07	0.000263	0.009674	0.008191
MF	GO: 0044389	Ubiquitin-like protein ligase binding	6	0.12	0.000207	0.038328	0.031622

Figure 3

Bubble chart of upregulated DEGs.

Figure 4

Bubble chart of downregulated DEGs.

Figure 5

Chord diagram of upregulated DEGs.

Figure 6

Chord diagram of downregulated DEGs.

Reactome pathway analysis

The upregulated DEGs were chiefly enriched in TP53, which regulates transcription of genes involved in G2 cell cycle arrest, condensation of prometaphase chromosomes, nuclear pore complex (NPC) disassembly, resolution of sister chromatid cohesion, and Golgi cisternae pericentriolar stack reorganization, and the downregulated DEGs were predominantly enriched in diseases associated with glycosaminoglycan metabolism, extracellular matrix organization, keratan sulfate degradation, diseases of glycosylation, and keratan sulfate biosynthesis. Table 3 lists the top 10 reactome pathways of the DEGs.

Table 3

The top 10 significantly enriched reactome pathways of the upregulated and downregulated DEGs.

Reactome ID	Description	Count	P value	FDR	Gene symbol
Downregulated
R-HSA-3560782	Diseases associated with glycosaminoglycan metabolism	8	6.07E-11	1.10E-08	FMOD, CSPG4, ACAN, LUM, PRELP, GPC1, OGN, BCAN
R-HSA-1474244	Extracellular matrix organization	15	9.43E-11	1.10E-08	FMOD, VCAM1, TIMP2, LTBP3, ITGB2, COL8A2, ACAN, FBLN1, LUM, PLEC, LAMB2, COL12A1, BCAN, MMP2, LTBP2
R-HSA-2022857	Keratan sulfate degradation	5	6.58E-09	5.10E-07	FMOD, ACAN, LUM, PRELP, OGN
R-HSA-3781865	Diseases of glycosylation	9	1.08E-07	6.29E-06	FMOD, CSPG4, THBS2, ACAN, LUM, PRELP, GPC1, OGN, BCAN
R-HSA-2022854	Keratan sulfate biosynthesis	5	4.70E-07	1.87E-05	FMOD, ACAN, LUM, PRELP, OGN
R-HSA-1630316	Glycosaminoglycan metabolism	8	4.83E-07	1.87E-05	FMOD, CSPG4, ACAN, LUM, PRELP, GPC1, OGN, BCAN
R-HSA-1638074	Keratan sulfate/keratin metabolism	5	1.30E-06	4.31E-05	FMOD, ACAN, LUM, PRELP, OGN
R-HSA-71387	Metabolism of carbohydrates	10	6.05E-06	0.000176	FMOD, CSPG4, ACAN, LUM, PRELP, GPC1, OGN, BCAN, PYGB, MAN2B1
R-HSA-3000178	ECM proteoglycans	5	7.15E-05	0.001849	FMOD, ACAN, LUM, LAMB2, BCAN
R-HSA-1474228	Degradation of the extracellular matrix	6	0.000144	0.003091	TIMP2, COL8A2, ACAN, COL12A1, BCAN, MMP2
Upregulated
R-HSA-6804114	TP53 regulates transcription of genes involved in G2 cell cycle arrest	4	1.37E-06	0.000294	AURKA, CDK1, BAX, CCNB1
R-HSA-2514853	Condensation of prometaphase chromosomes	3	1.68E-05	0.001539	CDK1, CCNB2, CCNB1
R-HSA-3301854	Nuclear pore complex (NPC) disassembly	4	2.21E-05	0.001539	CDK1, CCNB2, CCNB1, NEK7
R-HSA-2500257	Resolution of sister chromatid cohesion	6	2.88E-05	0.001539	CDK1, CENPM, CCNB2, CENPN, CCNB1, BIRC5
R-HSA-162658	Golgi cisternae pericentriolar stack reorganization	3	3.66E-05	0.001564	CDK1, CCNB2, CCNB1
R-HSA-202733	Cell surface interactions at the vascular wall	6	4.62E-05	0.001645	PSG7, PSG3, PSG8, PSG11, PSG4, PSG9
R-HSA-6791312	TP53 regulates transcription of cell cycle genes	4	8.53E-05	0.002274	AURKA, CDK1, BAX, CCNB1
R-HSA-68886	M phase	9	9.46E-05	0.002274	CDK1, CENPM, KIF20A, CCNB2, CENPN, UBE2C, CCNB1, BIRC5, NEK7
R-HSA-176412	Phosphorylation of the APC/C	3	9.58E-05	0.002274	CDK1, UBE2C, CCNB1
R-HSA-2980766	Nuclear envelope breakdown	4	0.000108	0.002303	CDK1, CCNB2, CCNB1, NEK7

The PPI network was screened, comprising 124 nodes and 355 edges (Figure 7). The top 10 hub genes were identified by cytoHubba plugin using the Maximal Clique Centrality (MCC) method, including TYMS, AURKA, CENPN, NUSAP1, CENPM, TPX2, CDK1, UBE2C, BIRC5, and CCNB1 (Table 4, Figure 8). Among these genes, CDK1 and CCNB1 showed the highest node degree. Then, 2 eligible modules of the gene interaction network were screened using the Molecular Complex Detection plugin (Figure 9). Module 1, consisting of 14 nodes and 91 edges, had 14 scores, while module 2, including 4 nodes and 6 edges, had 4 scores. Reactome pathways analysis revealed that genes in module 1 were primarily associated with M Phase, Resolution of Sister Chromatid Cohesion, Cell Cycle Checkpoints, Mitotic Prometaphase, and Condensation of Prometaphase Chromosomes, while in module 2 they were keratan sulfate degradation, keratan sulfate biosynthesis, keratan sulfate/keratin metabolism, diseases associated with glycosaminoglycan metabolism, and extracellular matrix organization. The top 10 reactome pathways of these modules are listed in Table 5.

Figure 7

PPI network of DEGs.

Table 4

The top 10 hub genes.

Gene symbol	Degree	MCC
CDK1	21	6.23E+09
CCNB1	21	6.23E+09
AURKA	18	6.23E+09
TPX2	15	6.23E+09
UBE2C	15	6.23E+09
CENPM	14	6.23E+09
CENPN	14	6.23E+09
NUSAP1	14	6.23E+09
BIRC5	14	6.23E+09
TYMS	14	6.23E+09

Figure 8

The top 10 hub genes.

Figure 9

(A, B) Two modules.

Table 5

The top 10 significantly enriched reactome pathways of module 1 and 2.

Reactome ID	Description	Count	P value	FDR	Gene symbol
Module 1
R-HSA-68886	M Phase	8	1.48E-09	4.09E-08	CDK1, UBE2C, CCNB1, CCNB2, BIRC5, KIF20A, CENPN, CENPM
R-HSA-2500257	Resolution of sister chromatid cohesion	6	2.16E-09	4.09E-08	CDK1, CCNB1, CCNB2, BIRC5, CENPN, CENPM
R-HSA-69620	Cell cycle checkpoints	7	7.98E-09	1.01E-07	CDK1, UBE2C, CCNB1, CCNB2, BIRC5, CENPN, CENPM
R-HSA-68877	Mitotic prometaphase	6	3.48E-08	3.30E-07	CDK1, CCNB1, CCNB2, BIRC5, CENPN, CENPM
R-HSA-2514853	Condensation of prometaphase chromosomes	3	1.81E-07	1.37E-06	CDK1, CCNB1, CCNB2
R-HSA-162658	Golgi cisternae pericentriolar stack reorganization	3	3.99E-07	2.52E-06	CDK1, CCNB1, CCNB2
R-HSA-6804114	TP53 regulates transcription of genes involved in G2 cell cycle arrest	3	8.91E-07	4.82E-06	CDK1, CCNB1, AURKA
R-HSA-176412	Phosphorylation of the APC/C	3	1.06E-06	5.01E-06	CDK1, UBE2C, CCNB1
R-HSA-69275	G2/M transition	5	1.45E-06	5.61E-06	CDK1, CCNB1, CCNB2, AURKA, TPX2
R-HSA-453274	Mitotic G2-G2/M phases	5	1.53E-06	5.61E-06	CDK1, CCNB1, CCNB2, AURKA, TPX2
Module 2
R-HSA-2022857	Keratan sulfate degradation	3	5.73E-09	1.81E-08	ACAN, FMOD, LUM
R-HSA-2022854	Keratan sulfate biosynthesis	3	6.56E-08	1.04E-07	ACAN, FMOD, LUM
R-HSA-1638074	Keratan sulfate/keratin metabolism	3	1.20E-07	1.26E-07	ACAN, FMOD, LUM
R-HSA-3560782	Diseases associated with glycosaminoglycan metabolism	3	2.13E-07	1.68E-07	ACAN, FMOD, LUM
R-HSA-1474244	Extracellular matrix organization	4	6.33E-07	4.00E-07	ACAN, FMOD, MMP2, LUM
R-HSA-3000178	ECM proteoglycans	3	1.40E-06	7.38E-07	ACAN, FMOD, LUM
R-HSA-1630316	Glycosaminoglycan metabolism	3	6.16E-06	2.78E-06	ACAN, FMOD, LUM
R-HSA-3781865	Diseases of glycosylation	3	9.47E-06	3.74E-06	ACAN, FMOD, LUM
R-HSA-71387	Metabolism of carbohydrates	3	8.32E-05	2.92E-05	ACAN, FMOD, LUM
R-HSA-1474228	Degradation of the extracellular matrix	2	0.001018	0.0003214	ACAN, MMP2

Discussion

OA is a disease involving the whole joint, and DMLs plays a critical role in the process of knee OA [7-9,29-31]. However, the treatment of DMLs remains a conundrum. Meniscectomy has no benefit over non-operative treatments, while non-operative measures such as nonsteroidal anti-inflammatory drugs, exercise therapy, and rest are the mainstay of management for DMLs [32]. Therefore, better understanding of the underlying pathogenesis of DMLs is of pivotal importance for identification of therapeutic targets and may drive future studies to search for biological treatments of this disease. The high-throughput microarray technology combined with bioinformatics analysis has been widely used in various diseases to predict potential molecular therapeutic targets. In the present study, the microarray dataset GSE52042, comprising healthier and diseased meniscus samples, was obtained, and a bioinformatics analysis was completed. A total of 154 DEGs were identified, including 64 upregulated and 90 downregulated DEGs. Furthermore, in order to probe into the biological significances behind these DEGs, GO, Reactome pathway, and PPI analysis were performed.

In view of the results of GO terms enrichment analysis, we linked the downregulated DEGs with ECM organization, collagen-containing ECM, and ECM structural constituent, and observed the upregulated DEGs had a close relationship with positive regulation of mitotic cell cycle, spindle, and ubiquitin-like protein ligase binding. The meniscus ECM, which plays a crucial role in sustaining structural integrity and mechanical properties of meniscus, is primarily composed of water, collagen, proteoglycan, glycoprotein, elastin, and linking protein [33,34]. Histopathological analysis showed that degenerated menisci could result in a loss of ECM structure [35]. Herwig et al. demonstrated that the collagen and glycosaminoglycan contents (chondroitin 4-sulphate, dermatan sulphate, and keratan sulphate) declined with increasing meniscus degeneration [36]. Meniscus cells are composed by different cell types, which possess proliferation ability and hasten meniscus repair [37-41]. Positive regulation of mitotic cell cycle can result in cell proliferation [42]. Nishida et al. observed fibroblasts and fibrochondrocyte-like cells proliferated but glycosaminoglycans and the cell count decreased remarkably after experimental dog meniscus tearing, indicating that ECM disorder and cell proliferation could coexist in the degenerative meniscus tear [43]. Lopez-Franco et al. showed that decreased proteoglycan, cartilage oligomeric matrix protein, and cell population altered the ECM organization in degenerative menisci, and actively dividing cells were found in all recently injured menisci samples and in 1 degenerative menisci sample [44].

The downregulated DEGs reactome pathways were predominantly involved in diseases associated with glycosaminoglycan metabolism, extracellular matrix organization, keratan sulfate degradation, diseases of glycosylation, and keratan sulfate biosynthesis. These results indicated the disorder of ECM was present in DMLs. Reactome pathway analysis of the upregulated DEGs were predominantly enriched in TP53, which regulates transcription of genes involved in G2 cell cycle arrest, condensation of prometaphase chromosomes, nuclear pore complex (NPC) disassembly, resolution of sister chromatid cohesion, and Golgi cisternae pericentriolar stack reorganization, and these were related to cell proliferation, consistent with the results of GO terms enrichment analysis of upregulated DEGs.

The PPI network was screened to predict the connections of protein functions of DEGs, and the top 10 hub genes – TYMS, AURKA, CENPN, NUSAP1, CENPM, TPX2, CDK1, UBE2C, BIRC5, and CCNB1 – were presented. Some hub genes are closely associated with osteoarthritic chondrocytes. Cyclin-dependent kinase 1 (CDK1), a highly conserved small protein, is a member of the cyclin-dependent kinase family that serves as a serine/threonine kinase and the catalytic subunit of mitosis-phase promoting factor, playing a pivotal role in controlling cell cycle, such as mitotic onset and driving G2-M transition [45-49]. CDK1 was detected in human and rat osteoarthritic chondrocytes [50]. Saito et al. demonstrated that overexpression of CDK1 could cause chondrocyte proliferation, and its low expression could result in chondrocyte differentiation and maturity, indicating CDK1 was essential for chondrocytes proliferation and differentiation [51]. Cyclin B1 (CCNB1) is one of cell cycle proteins which are essential regulators of cell division, and complexes with CDK1 to form a mitosis-phase-promoting factor that promotes mitosis [46,47]. Zhou et al. found that platelets could upregulate the proliferation of rat osteoarthritic chondrocytes via stimulation of the ERK/CDK1/CCNB1 pathway [52]. Baculoviral inhibitor of apoptosis repeat-containing 5 (BIRC5), also called survivin, is a protein that belongs to the inhibitor of apoptosis family, which localizes to the spindle, expressed in the G2/M phase, and plays a key role in the regulation of mitosis and protection of cells from apoptosis [53,54]. BIRC5 was detected in human osteoarthritic chondrocytes, and its suppression expression can lead to reduced rates of chondrocyte proliferation and G2/M blockade [55]. Aurora kinase A (AURKA) is a type of Aurora kinase and is crucial for the successful execution of mitosis [56]. AURKA was highly expressed in human osteoarthritic chondrocytes and overexpression of AURKA could contributes to the development of knee osteoarthritis [57]. In the present study, AURKA, CDK1, BIRC5, and CCNB1 were upregulated in diseased menisci compared to healthier menisci in late-stage knee OA, indicating diseased meniscus cells tended to cell proliferation.

Two modules were extracted from the PPI network. Most of the genes presented in module 1 are the aforementioned hub genes. The genes in module 2 relate to ECM metabolism. Matrix metalloproteinase-2 (MMP2) is a matrix metalloproteinase with ECM-degrading components [58]. Stone observed that late-stage OA meniscus cell cultures secreted more MMP2 than normal meniscus cell cultures when stimulated by pro-inflammatory factors [59]. Cook et al. found that when stimulated by interleukin-1β, MMP2 expression was decreased, but other matrix metalloproteinases were increased in canine meniscus cell cultures [60]. James et al. found that co-culture of meniscus and the infrapatellar fat pad, which is a major sources of pro-inflammatory factors, increased sulfated glycosaminoglycan degeneration but decreased MMP2 production [61]. In the present study, MMP2 was downregulated in diseased menisci compared to healthier menisci. Lumican (LUM) and fibromodulin (FMOD) are Class II small leucine-rich proteoglycans and are components of tendon and cartilage, with effects on collagen fibril structure and tissue stiffness [62,63]. Önnerfjord et al. found that LUM and FMOD were highly expressed in menisci, and the relative content between them in different tissues was negatively correlated [64]. Melrose et al. found that the expression of LUM and FMOD in menisci were consistent in late-stage knees OA versus age-matched normal knee joints [65]. Endo et al. observed that LUM was upregulated in menisci in early-phase knee OA [66]. Hellio observed that the FMOD expression was downregulated in the lateral meniscus at 3 weeks after anterior cruciate ligament transection for experimental OA in rabbits, but returned to normal at 8 weeks after surgery [67]. Aggrecan (ACAN) is a macromolecular proteoglycan and a crucial component of meniscus ECM, helping to maintain biomechanical properties of ECM [68-71]. Favero et al. found that OA meniscus and synovium could release inflammatory molecules to induce meniscus ACAN degradation and ECM loss [72]. Freymann et al. observed that human serum could stimulate meniscal cells obtained from advanced degenerative menisci in end-stage OA patients to express aggrecan [73]. In the present study, LUM, FMOD, and ACAN were downregulated in diseased menisci compared to healthier menisci in late-stage knee OA, indicating the ECM structure was impaired in diseased menisci.

According to the present results, some of the genes promoting cell proliferation were upregulated, while some of the genes constituting meniscus ECM were downregulated in the diseased meniscus samples, indicating that cell repair and ECM degradation may coexist in the process of osteoarthritic DMLs, which reflected the disorder of cell proliferation and the impairment of ECM integrality in osteoarthritic DMLs. Better understanding of these unique biological features may help develop potential therapies for osteoarthritic DMLs.

Conclusions

We conducted a comprehensive bioinformatics analysis on gene expression profile of OA meniscus. Key genes and pathways were identified to provide more detailed molecular mechanisms for the process of osteoarthritic DMLs and shed light on potential therapeutic targets. Nevertheless, further experiments are needed to further validate the identified genes and pathways.