
==== Front
iScience
iScience
iScience
2589-0042
Elsevier

S2589-0042(24)01966-7
10.1016/j.isci.2024.110741
110741
Article
The role of mitochondrial autophagy in osteoarthritis
Wang Genchun 123
Zhang Xiong 13
Xu Jingting 1
Hou Liangcai 1
Guo Zhou 1
Sun Kai 1085844308@qq.com
1∗
Guo Fengjing guofjdoc@163.com
14∗∗
1 Department of Orthopedics, Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, Hubei 430030, China
2 Orthopedic Medical Center, Union Hospital, Fujian Medical University, Fuzhou, Fujian 350000, China
∗ Corresponding author 1085844308@qq.com
∗∗ Corresponding author guofjdoc@163.com
3 These authors contributed equally

4 Lead contact

15 8 2024
20 9 2024
15 8 2024
27 9 11074118 11 2023
5 2 2024
12 8 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Summary

Osteoarthritis (OA) is a progressive degenerative joint disease, and the underlying molecular mechanisms of OA remain poorly understood. This study aimed to elucidate the relationship between mitochondrial autophagy and OA by identifying key regulatory genes and their biological functions. Utilizing bioinformatics analyses of RNA expression profiles from the GSE55235 dataset, we identified 2,136 differentially expressed genes, leading to the discovery of hub genes associated with mitochondrial autophagy and OA. Gene set enrichment analysis (GSEA) revealed their involvement in critical pathways, highlighting their potential roles in OA pathogenesis. Furthermore, our study explored the immunological landscape of OA, identifying distinct immune cell infiltration patterns that contribute to the disease’s inflammatory profile. We also evaluated the therapeutic potential of drugs targeting these hub genes, suggesting potential approaches for OA treatment. Collectively, this study advances our knowledge of mitochondrial autophagy in OA and proposes promising biomarkers and therapeutic targets.

Graphical abstract

Highlights

• Analyzed mechanisms of mitochondrial autophagy-related genes and osteoarthritis

• GSEA, WGCNA, GO, and KEGG analysis revealed biological processes

• BNIP3 potential therapeutic target for osteoarthritis

• Developing hub gene diagnostic model and immune correlation

Physiology; Cell biology; Bioinformatics; Transcriptomics

Subject areas

Physiology
Cell biology
Bioinformatics
Transcriptomics
Published: August 15, 2024
==== Body
pmcIntroduction

Osteoarthritis (OA), the most common degenerative joint disease globally, substantially impacts life quality with symptoms such as joint deformity and functional loss.1,2,3,4 Pathologically, OA is characterized by articular cartilage destruction, synovial hyperplasia, osteophytes, and subchondral bone sclerosis, yet its pathogenesis and molecular intricacies remain elusive.5 Current OA management primarily revolves around pharmacological interventions and surgery.6

Autophagy, a cellular homeostasis mechanism, plays a dual regulatory role in the progression of OA. It protects chondrocytes from oxidative stress in the early stages, while in the advanced stages, diminished autophagy contributes to cartilage degradation.7,8 Mitophagy, a subtype of autophagy crucial for chondrocyte mitochondrial health, involves key molecules such as PINK1, PRKN, BNIP3, and MFN2. The PINK1-PRKN pathway, particularly studied in nucleus pulposus cells (NP cells) and chondrocytes, offers promising therapeutic insights for OA,9 though its detailed regulatory mechanisms in OA require further exploration.

The advent of bioinformatics has greatly advanced protein research, aiding in the elucidation of protein functions, molecular pathways, and drug development.10 This approach has been pivotal in uncovering OA pathogenesis and identifying potential therapeutic targets. Key differentially expressed genes (DEGs) and biomarkers such as circRNAs and ferroptosis-related genes have been identified in OA through bioinformatics analyses,11,12,13 contributing to a growing repository of potential disease markers and potential therapeutic avenues.14,15,16

However, challenges like sample heterogeneity and confounding factors limit the accuracy of differential biomarkers. To mitigate these issues, a combination of DEGs, weighted co-expression networks (WGCNA), and validation across multiple independent datasets is recommended for comprehensive analyses. Our study utilizes microarrays and bioinformatics to construct a gene regulatory network around mitochondrial autophagy in OA.

Results

The identification of hub genes that regulate mitochondrial autophagy and OA

To investigate the relationship between mitophagy and osteoarthritis (OA), we conducted a comprehensive analysis (Figure 1). First, we identified 2,136 DEGs from GSE55235 using the screening conditions “p < 0.05 and |logFC|>1.5” (Table S1). A volcano plot shows these DEGs (Figure 2A). Heat maps for the top 20 DEGs were plotted (Figure 2B). We constructed a weighted gene co-expression network from GSE55235 after removing aberrant samples and filtering genes. Using a soft threshold force of 5, the scalar independence reached 0.87 and the average connectivity reached 69.86 (Figure 3A). As a result of dynamic tree-cutting with a sensitivity set to 3, a minimum module size of 30, and a module merging threshold of 0.25, 31 different co-expression modules were identified (Table 1) (Figure 3B), along with module feature vector clustering (Figure 3C). A correlation analysis was subsequently conducted between each module and clinical traits. There was the highest positive correlation between the MElightcyan1 module and OA (r = 0.97), while there was the highest negative correlation between the MEskyblue module and OA (r = −0.63) (Figure 3D). Here, the MElightcyan1 module, which contains 1,594 genes, has the largest absolute correlation and is selected for additional analysis to enrich the MElightcyan1 module genes. A total of 269 genes with high connectivity in the clinically significant module were identified as hub genes based on the cut-off criteria (| Module Membership [MM] | > 0.8). The Pearson correlation analysis revealed significant associations between the identified hub genes and OA, with p values less than 0.05. This statistical robustness underscores the strength of these genes’ correlations with OA pathogenesis. According to MM and Gene Signifificance (GS) correlation analysis, these genes may have a strong association with modules and phenotypes (COR = 0.98, p < 0.01) (Figure 3E). When the sensitivity is set to 3, the minimum module size is set to 30, and the module merging threshold is set to 0.25, 29 separate co-expression modules can be obtained by dynamic tree cutting, and a gene clustering dendrogram was drawn being dependent on the results (Figure 3F). In addition, we clustered the samples of the GSE55235 dataset to obtain a clustering dendrogram of the NC and OA groups (Figure 3G) Through the intersection of 2,136 differential genes (Tables S1 and S4), 269 MElightcyan1 module hub genes (Table S2), and 2,364 mitochondrial autophagy-related genes screened (Table S3) by GSE55235, 31 hub genes related to mitochondrial autophagy and OA was obtained (Figure 4A). The STRING database was used to construct a Protein-Protein Interaction Networks (PPI) network that was used to reveal the interactions of each protein. The PPI network was actuated by the degree of target connectivity (DEgree) from largest to smallest (Figure 4B). Four cytoHubba algorithms (degree, Maximum Neighborhood Component [MNC], stress, and Maximal Clique Centrality [MCC]) were consistently able to determine hub genes. The six genes with the highest score were thought hub genes for OA (Table 2), namely MYC, BNIP3, GABARAPL1, PPARGC1A, MCL1, and BAG3 (Figures 4C and 4D). The violin plot shows that in GSE55235, the expression levels of MYC, BNIP3, GABARAPL1, MCL1, and BAG3 are significantly decreased in the OA group compared to the normal group (p < 0.05) (Figure 5A). Confirming the relevance of these observed results highlights their potential roles in OA. Similarly, we evaluated the mRNA expression levels of these six genes in the GSE55457 dataset and obtained partially similar results to those mentioned previously (Figure S1).Figure 1 Flowchart of the study

Figure 2 Genes differentially expressed between NC and OA samples

(A) Identification of differentially expressed genes. The volcano plot shows the differentially expressed genes in the GSE55235 dataset, green indicates genes highly expressed in NC, red indicates genes highly expressed in OA, and gray indicates genes with no significant change.

(B) The heatmap displays the top twenty upregulated and downregulated differentially expressed genes in NC or OA samples.

Figure 3 Results of WGCNA

(A) Fit indices of the corresponding scale-free topology models under different soft thresholds.

(B) The corresponding average connectivity under different soft thresholds.

(C) Eigenvector clustering of different modules.

(D) Heatmap of the correlation of different modules with phenotypes.

(E) Correlation of module membership and gene significance in the lightcyan1 module.

(F) Clustering dendrogram of genes.

(G) Clustered dendrograms of different samples.

Table 1 Identification of co-expression modules

Module name	Number	
bisque4	98	
blue	870	
brown	727	
brown4	33	
cyan	204	
darkgreen	75	
darkgrey	68	
dark magenta	126	
darkolivegreen	60	
darkorange2	34	
darkslateblue	30	
darkturquoise	71	
floralwhite	36	
ivory	838	
lightcyan1	1594	
lightsteelblue1	40	
lightyellow	178	
mediumpurple3	40	
orangered4	631	
paleturquoise	61	
plum1	46	
saddlebrown	62	
sienna3	55	
skyblue	62	
skyblue3	50	
steelblue	61	
tan	147	
violet	60	
yellowgreen	55	

Figure 4 Acquisition and correlation analysis of hub genes

(A) The intersection of 31 hub genes was obtained through the analysis of DEGs, WGCNA of the lightcyan1 module genes, and mitochondrial autophagy-related genes.

(B) Ranking of the connections of 31 central genes. The central genes with lower rankings are hidden, and based on the rankings, the 6 hub genes with the most closely related connections were selected.

(C) The identification of six central genes using four algorithms in CytoHubba and the evaluation of their protein-protein interaction patterns were conducted.

(D) Interaction of 6 hub genes obtained by string database.

Table 2 The differential expression results of 6 single genes in the NC and OA groups

Label	Comparison group (mean ± SD)	control group (mean ± SD)	Mann-Whitney test	
MYC	NC(10.48 ± 0.59)	OA(7.95 ± 0.57)	1.1e-5	
BNIP3	NC(8.06 ± 0.58)	OA(6.45 ± 0.50)	1.1e-5	
GABARAPL1	NC(8.93 ± 0.61)	OA(7.57 ± 0.42)	7.6e-5	
PPARGC1A	NC(4.07 ± 1.05)	OA(3.86 ± 0.83)	0.58	
MCL1	NC(9.79 ± 0.37)	OA(8.51 ± 0.28)	1.1e-5	
BAG3	NC(10.11 ± 0.74)	OA(8.20 ± 0.27)	1.1e-5	

Figure 5 Expression of hub genes in articular cartilage of NC and OA samples in the GSE55235 dataset

Biological processes and pathways of hub gene enrichment

Using enrichment analysis, we were able to understand the potential biological functions of these genes. Gene ontology (GO) analysis included the analysis of BP, CC, and MF. The results showed that these six genes are mainly engaged in autophagy, the process utilizing autophagic mechanism (Figure 6A), neuron projection, mitochondrion (Figure 6B), transcription factor binding, and Tubin binding (Figure 6C). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed that five out of six genes are involved in autophagy-related biological processes, including mitophagy, autophagy, FoxO, and Apelin signaling pathways (Figure 6D). The gene set enrichment analysis (GSEA) results indicate that these six genes were related to the p53 signaling pathway, NOD-like receptor signaling pathway, lysosome, glycerol lipid metabolism, antigen processing, and presentation (Table 3). GSEA indicated significant associations of these genes with the p53 signaling and NOD-like receptor signaling pathways, with Normalied Enrichment Score (NES) values >1.5, Nominal (NOM) p values <0.05, and false discovery rate (FDR) q values <0.25, suggesting their roles in OA progression (Figure 6E).Figure 6 Functional correlation analysis

Hub genes were analyzed by GO, KEGG, and GSEA.

(A) Biological process (BP).

(B) Cellular component (CC).

(C) Molecular function (MF).

(D) KEGG signaling pathway.

(E) GSEA shows the enriched signaling pathways associated with the hub genes that were identified through the filtering. ES, enrichment score; NP, Nominal p value, ∗P < 0.05, ∗∗p < 0.01.

Table 3 GSEA analysis

Term	ES	NES	p value	FDR	FWER	
LYSOSOME	−0.5120	1.7271	0.0482	0.5739	0.2610	
GLYCEROLIPID_METABOLISM	−0.5465	−1.4643	0.0261	0.6429	0.9020	
ANTIGEN_PROCESSING_AND_PRESENTATION	−0.4994	−1.4782	0.0380	0.6595	0.8860	
NOD_LIKE_RECEPTOR_SIGNALING_PATHWAY	0.4932	1.5230	0.0103	0.6892	0.8090	
P53_SIGNALING_PATHWAY	0.3910	1.3723	0.0467	1.0000	0.9680	

RT-qPCR validation of the hub genes

RT-qPCR experiments were carried out to verify the bioinformatics results. The results showed that the expression of MYC, BNIP3, GABARAPL1, BAG3, MCL1, and PPARGC1A were markedly different in the normal and IL-1β-induced groups, indicating that the bioinformatics data mining results are more reliable. RT-qPCR validation demonstrated significant expression differences for hub genes between normal and OA groups, with p values <0.05 and fold changes >2. The experimental design included triplicate measurements from samples per group, enhancing the validation’s reliability (Figure 7A).Figure 7 Validation of Hub genes at the cytological and histological levels

(A) RT-PCR results showed mRNA expression levels of Hub genes in mouse chondrocytes in the NC group versus the IL-1β group.

(B) Protein expression of BNIP3 in the NC group versus the TBHP-intervened chondrocyte group.

(C) Immunofluorescence of BNIP3 protein expression in the NC group versus the TBHP-intervened chondrocyte group (bar: 50 μm).

(D) The staining of guanine green and OARSI scores in mouse knee cartilage (bar: 50 μm).

(E) Immunohistochemical staining of BNIP3 in mouse cartilage (bar: 50 μm).

All data represent mean ± SD. ∗p < 0.05, ∗∗p < 0.01.

Decreased expression of BNIP3 in osteoarthritis

BNIP3 was selected for further analysis for two reasons. As BNIP3 is centrally located in the hub gene network, it regulates other genes; and the HIF-1α-BNIP3 axis is involved in mitochondrial autophagy. It mitigates OA progression.17 Western blot analysis quantitatively showed a marked decrease in BNIP3 protein levels in OA samples compared to normal, with a significant statistical difference (p < 0.05), supporting BNIP3’s role in OA pathogenesis (Figure 7B). Immunofluorescence staining and quantitative analysis showed that the protein expression of BNIP3 in chondrocytes was down-regulated in the inflammatory environment of OA, which is consistent with the decreased mRNA expression of BNIP3 in RT-qPCR (p < 0.05) (Figure 7C). The differential expression of BNIP3 protein in the normal group and OA group may be attributed to TBHP promoting the enhancement of chondrocyte oxidative stress capacity, leading to increased Reactive oxygen species (ROS) production, more severe mitochondrial damage, and rapid onset of mitochondrial dysfunction and inhibition of BNIP3 expression.

We assessed the severity of OA using the Osteoarthritis Research Society International (OARSI) score and rated cartilage changes using Safranin-O solid green staining in a surgically induced destabilizing the medial meniscus (DMM) mouse model. A sizable amount of cartilage degradation was observed in the DMM group. OA severity, assessed using OARSI scores, showed significant cartilage degradation in the DMM group (p < 0.01). Safranin-O staining quantitatively corroborated these findings (Figure 7D). In the routine and DMM groups, immunohistochemical staining analysis revealed a decreased ratio of BNIP3-positive cells in the articular cartilage of DMM mice compared to the control group, which was consistent with the results of in vitro experiments. These results demonstrate the presence of reduced BNIP3 expression after articular cartilage damage (p < 0.05) (Figure 7E).

Developing and validating diagnostic models

A logistic regression algorithm based on GSE55235 was utilized to construct a multi-gene prediction model. The prognostic significance of these genes MYC, BNIP3, GABARAPL1, PPARGC1A, MCL1, and BAG3 in 20 samples was assessed using the cox method (Table 4). The overall prognosis was significantly changed (log-test = 0.0009, sc-test = 0.002, wald-test = 0.493, C-index = 1). The logistic regression model, based on GSE55235, exhibited excellent predictive capability for OA with an Area Under Curve (AUC) value of 0.95. Specific logistic regression parameters were optimized using stepwise regression (Figure 8A). The results showed that our constructed osteoarthritis prediction model had a good predictive effect. Next, we used GSE55457, GSE82107, and GSE169077 as validation datasets for further validation. The results showed that the AUC = 0.99 for the GSE55457 model (Figure 8B) and AUC = 0.94 for the integrated dataset model of GSE82107 and GSE169077 (Figure 8C), indicating that the model has some guiding significance for diagnosing OA patients in clinical applications. Validation using GSE55457, GSE82107, and GSE169077 demonstrated consistent predictive performance with AUC values ranging from 0.94 to 0.99, underscoring the model’s clinical relevance. Detailed dataset analysis including sample sizes and statistical significance of model performance is provided.Table 4 Muticox analysis

Tag	exp(coef)	Pr(>|z|)	lower 0.95	upper 0.95	coef	
MYC	0.791001953782197	0.73215495616388	0.206580043856656	3.02877315352595	0.234454841202	
BNIP3	0.841898652166623	0.852474951548842	0.13725446763787	5.16408210762249	−0.1720956375821	
GABARAPL1	0.358882992436383	0.378381420922189	0.0367047960585653	3.50899653698081	−1.0247588699804	
PPARGC1A	0.945637570779286	0.914300181129945	0.341672547508795	2.61721470393032	−0.0558959009080	
MCL1	0.784608184676227	0.87195619850546	0.0410747848987149	14.9875405307402	−0.2425708136064	
BAG3	0.223766912891137	0.274510210133301	0.0152605125092237	3.28112383346008	−1.4971503363435	

Figure 8 ROC curves and corresponding AUC values for the three diagnostic models

(A) Diagnostic models based on GSE55235.

(B) Diagnostic models based on GSE55457.

(C) Diagnostic models based on integrated datasets of GSE82107 and GSE169077.

Immune infiltration and immune-related factors

Clinical treatment sensitivity and disease diagnosis are a function of the immunological microenvironment. The CIBERSORT algorithm was used to evaluate the proportion of 22 immune cells in 10 NC samples and 10 OA samples. CIBERSORT analysis indicated significant differences in immune cell composition between OA and NC samples, providing detailed percentages of each immune cell type and their statistical differences (Figure 9A). Furthermore, we compared immune cell infiltration in the OA and NC samples in the boxplot, as well as the correlation between immune cells (Table 5). The results of immune cell infiltration showed that the percentage of B cells memory, plasma cells, T cells CD4 naive, and mast cells resting in the OA group was significantly higher than that of the NC group; while B cells naive, T cells CD4 memory resting, natural killer (NK) cells active, and the percentage of B cells naive, T cells CD4 memory resting, NK cells activated, and mast cells activated were significantly lower than the NC group (Figure 9B). Comparative analysis of immune cell types showed significant correlations among different immune cells in chondrocytes (p < 0.05) (Figure 9C). After further exploring the relationship between hub genes and immune infiltration, we conducted a correlation analysis of gene expression levels and immune cell contents of the selected hub genes. The results showed a strong association between hub genes and immune infiltration. Specifically, we found that MYC, BNIP3, GABARAPL1, MCL1, and BAG3 were highly positively correlated with resting CD4 memory T cells and enabled mast cells. Meanwhile, MYC, BNIP3, GABARAPL1, MCL1, and BAG3 were significantly and negatively correlated with plasma cells, B cells memory, and mast cells resting (Figure 9D).Figure 9 Analysis of immune infiltration

(A) The proportion of 22 immune cells corresponding to each sample in the dataset GSE55235.

(B) The proportion of immune cells in the NC and OA groups.

(C) Correlation of immune cells.

(D) Relationship between hub genes and immune cells.

Table 5 Immune components results

Label	Comparison group (mean ± SD)	Control group (mean ± SD)	Mann-Whitney test	
B cells naive	NC(0.03 ± 0.03)	OA(8.6e-3±0.01)	0.03	
B cells memory	NC(0.02 ± 0.03)	OA(0.07 ± 0.06)	0.02	
Plasma cells	NC(6.8e-4±1.8e-3)	OA(0.04 ± 0.05)	1.0e-3	
T cells CD8	NC(0.01 ± 0.02)	OA(0.02 ± 0.02)	0.28	
T cells CD4 naive	NC(0.0e+0 ± 0.0e+0)	OA(9.9e-3±0.02)	0.03	
T cells CD4 memory resting	NC(0.16 ± 0.05)	OA(0.02 ± 0.02)	1.8e-4	
T cells CD4 memory activated	NC(6.8e-3±0.01)	OA(0.01 ± 0.01)	0.07	
T cells follicular helper	NC(0.02 ± 0.01)	OA(0.01 ± 0.01)	0.38	
T cells regulatory (Tregs)	NC(0.01 ± 0.01)	OA(0.02 ± 0.01)	0.38	
T cells gamma delta	NC(0.01 ± 0.02)	OA(0.01 ± 0.01)	0.23	
NK cells resting	NC(3.4e-3±0.01)	OA(4.6e-3±9.4e-3)	0.39	
NK cells activated	NC(0.04 ± 0.04)	OA(6.1e-3±9.6e-3)	9.7e-3	
Monocytes	NC(0.04 ± 0.03)	OA(0.02 ± 0.03)	0.26	
Macrophages M0	NC(0.01 ± 0.02)	OA(0.08 ± 0.10)	0.19	
Macrophages M1	NC(0.06 ± 0.06)	OA(0.04 ± 0.02)	0.58	
Macrophages M2	NC(0.31 ± 0.12)	OA(0.39 ± 0.10)	0.14	
Dendritic cells resting	NC(0.03 ± 0.03)	OA(0.03 ± 0.02)	0.55	
Dendritic cells activated	NC(1.1e-3±3.4e-3)	OA(0.0e+0 ± 0.0e+0)	0.37	
Mast cells resting	NC(0.03 ± 0.05)	OA(0.17 ± 0.06)	3.7e-4	
Mast cells activated	NC(0.17 ± 0.12)	OA(2.0e-3±5.2e-3)	1.1e-4	
Eosinophils	NC(3.3e-3±7.6e-3)	OA(0.0e+0 ± 0.0e+0)	0.17	
Neutrophils	NC(0.01 ± 0.01)	OA(0.01 ± 0.01)	0.73	

Analysis of drug-gene interactions, related-proteins subcellular localization, and molecular dockings

A protein’s subcellular localization identifies its biological function. The protein subcellular localization of MYC, BNIP3, GABARAPL1, PPARGC1A, MCL1, and BAG3 was predicted using public compartments database. The subcellular localization of key proteins was predicted with high confidence scores, indicating their functional relevance in specific cellular compartments (Figure 10A). The development of potential therapeutic agents targeting genes could serve as a specific therapeutic strategy. Using data from DrugBank and Metascape databases,18 two candidate drugs were identified. Acetylsalicylic acid is a drug that targets MYC, and isosorbide targets MCL1 (Table 6). Acetylsalicylic acid is associated with stimulating autophagy as well as inhibiting apoptosis. Isosorbide is thought to inhibit cellular senescence and maintain normal oxidative stress homeostasis. Autodock Vina calculated the binding energies of the two candidate drugs with two proteins, as well as their binding poses and interactions. The binding energy is defined as the reduction in energy required for the entire system after the protein and molecule are bound. Generally, a binding energy of −4.25 kcal/mol indicates a certain interaction between the receptor protein and the ligand. The receptor (PDB structure of MYC is 6QB6, and the PDB structure of MCL1 is 6E16) was docked with the respective ligands. The protein and ligand files were converted to PDBQT format, water molecules were removed, and polar hydrogen atoms were added. The docking pocket was set to the entire protein, with a grid point spacing of 0.05 nm. The molecular binding sites between MYC, MCL1, and approved small molecule drugs were calculated separately, and the results showed a binding energy of −5.858 kcal/mol for MYC and acetylsalicylic acid (Figure 10B), and a binding energy of −4.445 kcal/mol for MCL1 and mannitol (Figure 10C), indicating a highly stable binding of the corresponding small molecule drugs with the targets.Figure 10 Protein subcellular localization and molecular docking structure of hub genes

(A) Protein subcellular localization of hub genes.

(B) Molecular docking structure of MYC and acetylsalicylic acid, hydrogen bonding of the binding sites is TYR-385.

(C) Molecular docking structure of MCL1 and isosorbide, hydrogen bonding of the binding sites is ARG-263.

Table 6 Molecular docking

Gene Symbol	Gene ID	Description	Subcellular location (Protein Atlas)	Drug (DrugBank)	
MYC	4609	MYC proto-oncogene, bHLH transcription factor	Nucleoplasm (enhanced)	Acetylsalicylic acid	
BNIP3	664	BCL2 interacting protein 3		
GABARAPL1	23710	GABA type A receptor-associated proteins like 1	Vesicles (approved)		
PPARGC1A	10891	PPARG coactivator 1 alpha	Nucleoplasm (supported)		
MCL1	4170	MCL1 apoptosis regulator, BCL2 family member	Mitochondria (enhanced)	Isosorbide	
BAG3	9531	BAG cochaperone 3	Cytosol (supported)		

Discussion

Osteoarthritis is a degenerative disease that results in the breakdown of articular cartilage, thickening of subchondral bone, formation of bone hyperplasia, synovitis, and joint collapse.19 The pathological process involves increased cartilage catabolism, decreased cartilage anabolism and repair, hypertrophy and death of chondrocytes, impaired autophagy function, osteoclast-mediated bone remodeling, and infiltrating and activating immune cells.9 Inflammation in OA differs from other autoimmune diseases as it is primarily mediated by the innate immune system. The interaction between the immune system and local tissue damage and metabolic dysfunction further complicates the disease process.20

Our study extends these findings by elaborating on the specific molecular pathways involved in OA, such as mitochondrial autophagy and its regulation by key genes. This approach bridges a critical gap in understanding the complex interplay between genetic factors and OA pathogenesis. Despite the expected increase in OA incidence with an aging population, the underlying pathogenesis remains poorly understood, and no disease-modifying drugs are currently available.21 Our study aimed to investigate the expression pattern and active process of mitochondrial autophagy-related genes in OA using comprehensive and effective bioinformatics methods, to identify new targets for the diagnosis, treatment, and prognosis of early- to mid-stage OA. While the molecular mechanism of mitochondrial autophagy and the immune environment in OA is not fully understood, our study offers new insights that may pave the way for future therapeutic interventions.

Our study aimed to investigate the relationship between mitochondrial autophagy and OA progression. By analyzing RNA expression profiles of OA, we identified 2,159 DEGs and screened 6 genes related to mitochondrial autophagy and OA. Through GO analysis, we found these genes were involved in mitochondrial autophagy, cell homeostasis regulation, and apoptosis. Selective autophagy has been considered an effective approach to eliminating damaged cells under oxidative stress.22 Previous studies have shown that mitochondrial dysfunction is associated with cartilage matrix breakdown and inhibition of cartilage matrix synthesis, which are correlated with OA.23 Mitophagy is a fundamental mechanism for maintaining mitochondrial homeostasis and can be used to clear damaged mitochondria, avoid inflammation activation, and inhibit apoptosis.24 Studies have demonstrated that trehalose promotes the increase of autophagic flux, which can alleviate mitochondrial dysfunction in chondrocytes and relieve the progression of OA. This process is linked with BNIP3-mediated mitotic activity.25 Moreover, the PPARGC1A gene is essential in regulating mitochondrial biogenesis, while REDD1 controls the mitochondrial biogenesis of articular chondrocytes by regulating PPARGC1A to alleviate OA progression.26 Studies pointed out that MYC, BNIP3, BAG3, MCL1, PPARGC1A, and GABARAPL1 are commonly involved in regulating mitochondrial function and delaying the progression of degenerative diseases. Abnormal mitochondrial autophagy may be the key link between autophagy function impairment and extracellular matrix (ECM) homeostasis disorder in several degenerative diseases.27,28,29,30 Regarding the clinical applications of our findings, the identified hub genes not only offer a deeper understanding of OA’s molecular underpinnings but also serve as promising biomarkers for early detection and therapeutic targets. The differential expression and pathway involvement of these genes suggest their potential role in diagnostic panels or as candidates for potential drug development aimed at modulating mitochondrial autophagy pathways to halt or reverse OA progression.

In contrast to rheumatoid arthritis (RA), synovial immune activation in OA exhibits distinctive characteristics, aiding in the differentiation between mild RA and inflammatory OA.31 This distinction is critical for patient stratification, diagnosis, and guiding treatment interventions.32 Early stages of OA are marked by the infiltration of inflammatory Th1/17 cells, underscoring a complex interplay between immune cells and humoral immunity, which is believed to drive the pathogenesis of autoimmune arthritis and related diseases.33 The cytokine profile in synovial fluid mirrors the effector and regulatory T cell composition within the synovium,34 with elevated expression of CCR5, CCR3, and CD161 among CD4+T cells in the joints indicating significant immune activation.35 Moreover, mast cells within the synovial tissue, known for their role in structural damage through the release of mediators such as tryptase, contribute to inflammation, cartilage degradation, and chondrocyte apoptosis in OA.36 Our investigation into immune infiltration disparities between normative and OA cohorts reveals significant variances in the prevalence of plasma cells, CD4+T cells, NK cells, and mast T cells in OA patients, aligning with findings from prior studies.37,38,39 This analysis illuminates the nuanced immune landscape of OA, offering insights into the disease’s immunopathology.

We validated the expression levels of hub genes in normal cartilage and OA models by RT-qPCR, western blot, immunofluorescence, and immunohistochemistry. Significant differences in mRNA expression levels of MYC, BNIP3, MCL1, GABARAPL1, BAG3, and PPARGC1A were observed between the normal and OA groups. Additionally, we address potential confounding factors and alternative explanations for the observed correlations, enhancing the scientific rigor of our study. Finally, we retrieved two drugs targeting the aforementioned genes from DrugBank and Metascape databases. Acetylsalicylic acid (aspirin) can inhibit the activity of cyclooxygenase (COX).40 Acetylsalicylic acid active metabolite, salicylate, stimulates autophagy flux and slows down aging by inhibiting acetyltransferase EP300,41 additionally, it reduces chondrocyte apoptosis and cartilage degradation in knee osteoarthritis by inhibiting MMPs, NO, and NF-κB signaling.42,43 Isosorbide can be used to correct the damaged redox balance between superoxide and nitric oxide, improve calcium circulation and contractility by reducing sarcoplasmic Ca2+ leakage, and alleviate vascular sclerosis caused by mitochondrial oxidative stress. It is also used for the treatment of congestive heart failure.44 Studies have shown that isosorbide can inhibit the senescence of mesenchymal stem cells, reduce senescence-related galactosidase (SA-β-gal) activity, and p21 expression.45 In addition, we conducted molecular docking experiments to confirm the binding sites of acetylsalicylic acid and nitrobenzene with their respective key proteins MYC and MCL1, with corresponding binding hydrogen bonds TYR-385, ARG-263. The potential drugs showed high binding energy with the key targets, stability, and potential to reduce off-target effects. These findings can serve as valuable resources for further drug development and clinical applications. Our molecular docking results provide detailed insights into the interaction between these drugs and target proteins, confirming their potential therapeutic applications in OA.

A limitation of this study is that gender differences were not fully controlled, which could result in some genes being more significantly expressed in females than males in the gene expression analysis. To mitigate this limitation and enhance the reliability of our findings, future studies will incorporate more rigorous control of gender variables, allowing for a more nuanced understanding of how sex differences may influence gene expression and OA progression.

These results not only highlight the clinical relevance of our findings but also delineate a clear path for future research, particularly in exploring the translational potential of these discoveries into clinical applications. By elucidating the complex molecular mechanisms underlying OA and proposing potential therapeutic targets, our study lays the groundwork for developing innovative treatment strategies. Our results suggest a pronounced association between mitochondrial autophagy and OA progression, revealing the significance of specific hub genes such as MYC, BNIP3, GABARAPL1, MCL1, PPARGC1A, and BAG3. These genes, as evidenced by their differential expression and involvement in key biological processes and pathways, may contribute to the development and exacerbation of OA through their roles in mitochondrial function and cellular stress responses.

Limitations of the study

However, this study also has some limitations. Firstly, the sample study did not exclude the gender factor, which may affect the universality of the results. Secondly, the data mainly came from articular cartilage tissue, lacking validation for surrounding tissues, which may limit the application of the diagnostic model and require further validation through clinical trials to determine the value of the diagnostic model. Finally, further validation of the expression of related genes can be conducted in other animal models.

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Professor Fengjing Guo (guofjdoc@163.com).

Materials availability

This study did not generate new unique reagents.

Data and code availability

• This study analyzes existing, publicly available data. The sources for the datasets are listed in the key resources table.

• This study does not report the original code. All codes were used in this study in alignment with recommendations made by authors of R packages in their respective user’s guide, which can be accessed at https://www.r-project.org/. All code used in the analyses is deposited on https://www.jianguoyun.com/p/DUTcSZsQjf26Cxjtr_oEIAA.

• Any additional information required to reanalyze the data used in this study is available from the lead contact upon request.

• Publicly available datasets were analyzed in this study. These data can be found Database: GSE55235, GSE55457, GSE82107, GSE169077.

Acknowledgments

This study was supported by the 10.13039/501100001809 National Natural Science Foundation of China (82372475 and 82172498 ). The PCR and western blot experiments were performed at the Experimental Medicine Research Center, Tongji Hospital, Huazhong University of Science and Technology. Dr. Genchun Wang and Dr. Xiong Zhang served as guarantors of this study, ensuring the integrity and accuracy of the data.

Author contributions

G.W. and X.Z. conceived and planned the study, and supervised the analysis work. G.W., X.Z., and J.X. conducted the analysis and drafted the initial manuscript. Z.G. and L.H. critically reviewed the integrity and credibility of the data analysis. K.S. and F.G. revised the manuscript and were responsible for the integrity of data acquisition and statistical analysis. All authors agreed to submit this manuscript, read and approve the final version, and take full responsibility for its content, including the accuracy of the data and the fidelity of the experiment to the registration agreement and its statistical analysis. Each author has unrestricted access to all data.

Declaration of interests

The authors declare no competing interests.

STAR★Methods

Key resources table

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
GAPDH	Proteintech Group, USA	60004-1-Ig; RRID: AB_2107436	
BNIP3	Cell Signaling Technology, USA	#44060	
Interleukin-1β (IL-1β)	R&D systems	#401-ML	
Tert-Butyl hydroperoxide solution (TBHP)	Sigma-Aldrich ,St Louis, MO, USA	418064	
Trypsin	Biosharp Life Sciences, Hefei, Anhui, China	BL527A	
Collagenase II	Biosharp Life Sciences, Hefei, Anhui, China	BS164	
4% Paraformaldehyde	Biosharp Life Sciences, Hefei, Anhui, China	BL539A	
DMEM/F12 medium	Hyclone, Logan, UT, USA	SH30023.01	
Fetal bovine serum Excellent	NEWZERUM LIMITED, Upper Riccarton, Christchurch, New Zealand	E500	
total RNA extraction kit	Toyobo, Japan	NPK-101	
First Strand cDNA Synthesis kit	Toyobo, Japan	RMD-101T	
SYBR Green real-time PCR Master Mix	Toyobo, Japan	QPK-201	
RIPA lysis buffer	Boster, Wuhan, Hubei, China	AR0102	
Protease inhibitors	Boster, Wuhan, Hubei, China	AR1182	
Phosphatase inhibitors	Boster, Wuhan, Hubei, China	AR1183	
BCA assay kit	Boster, China,	AR0146	
DAPI reagent	Boster, Wuhan, Hubei, China	AR1176	
TritonX-100	Biofroxx, Germany	143306	
DAB staining	Sigma-Aldrich ,St Louis, MO, USA	D8001-5G	
EDTA	EDTA decalcifying solution	E1171	
Safranin O Stain Solution	Solarbio Life Sciences, Beijing, China	G1067	
Fast Green FCF Stain Solution	Solarbio Life Sciences, Beijing, China	G1661	
PVDF membranes	Millipore, USA	IPVH00010	
SDS-PAGE gels	BioRad, Hercules, CA, USA	4561036	
	
Experimental models: Organisms/strains	
	
C57BL/6J mouse	GemPharmatech, Nanjing, China	Strain NO.N000013	
	
Oligonucleotides	
	
Primer for RT-qPCR, see Material and method	Tsingke Biotech,China	https://www.tsingke.com.cn/	
	
Software and algorithms	
	
ImageJ Version 1.52i	National Institutes of Health	https://www.nih.gov/	
R Software 4.1.1.	The R Project	https://www.r-project.org	
Pearson	Corrplot (The R package)	https://cran.r-project.org/web/packages/corrplot/index.html	
GraphPad Prism Version 7	GraphPad Software	https://www.graphpad.com/	
Adobe Photoshop CC Version 19.0	Adobe	https://www.adobe.com/	

Experimental model and study participant details

Animal experiments and OA model

The Laboratory Animal Center (Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China) approved all animal experimental procedures. (TJH-202201013). This study utilized twelve 7-week-old male C57BL/6J mice. In the sham-operated group, 6 mice were randomly selected. The sham-operated group underwent sham surgery after receiving an intrapulmonary injection of pentobarbital (35 mg/kg). The OA model was then established by destabilizing the medial meniscus (DMM) in 6 mice. The knee joints were gathered and fixed in 4% paraformaldehyde for further experiments after the mice were sacrificed 8 weeks later. After being stained histologically and immunohistochemically, knee joints were decalcified with 10% EDTA solution for two weeks and integrated into paraffin. We stained specimens with Safranin-O/fast green after cutting them in the sagittal plane to 5mm thicknesses. OARSI (Osteoarthritis Research Society International) scores were used to measure OA progression. In addition, further sections were deparaffinized, antigen extracted, and incubated with anti-BNIP3 antibody, followed by goat anti-rabbit secondary antibody. DAB staining and hematoxylin counterstaining were made on the sections. Under 200x magnification, images were taken with a fluorescence microscope (Evosflauto. Life Technologies, USA) and the number of BNIP3-positive cells was counted.

Method details

Microarrays Data Source and Differential Expression Analysis

We sourced four datasets (GSE55235,46 GSE55457,47 GSE82107,48 GSE169077,49) from the GEO database. Selected datasets comprised RNA expression profiles from OA and normal samples. Differential expression analysis was conducted using the R "limma" package with |log2FC| > 1.5 and p < 0.05 as cutoff values. Data visualization employed "pheatmap" and "ggplot2" packages for depicting DE genes (Table S4).

Mitochondrial autophagy-related genes and WGCNA

Mitochondrial autophagy-related genes (MRGs) were obtained from the Genecards Database,50 totaling 2364 genes after removing overlaps. WGCNA, performed using the R "WGCNA" package, analyzed the correlation between gene expression patterns and OA phenotypes.51 We applied dynamic tree-cutting methods for module detection and correlated these with OA traits using module Eigengenes.

Hub Gene Identification and Validation

The intersection of DE genes, WGCNA, and MRGs was examined using the R "VennDiagram" package to identify OA-related hub genes. Expression differences were visualized using violin plots, and the cytoHubba plugin estimated hub gene networks (Table S5). Statistical analyses, including t-test and Mann-Whitney U test, were employed to validate the findings.

Functional enrichment analysis and GSEA

Functional enrichment, encompassing GO and KEGG analysis, used org.Hs.Eg.Deb (R package) for gene annotations, and KEGG for pathway analysis (Table S6). ClusterProfiler package conducted enrichment analysis with p < 0.05 as the significance threshold. GSEA,52 based on pivotal gene expression, categorized genes by logFC and utilized the "enrichplot" package for result visualization.

Obtaining and culturing mouse chondrocytes

Previously described techniques were utilized to isolate and culture chondrocytes.53 Cartilage was harvested from five-day-old C57BL/6J mice. The primary chondrocytes were resuspended in DMEM/F12 medium, supplemented with 10% fetal bovine serum (FBS), 1% penicillin, and 1% streptomycin sulfate, and cultured at 37°C in an atmosphere containing 5% CO2. Experiments were conducted using the second-generation chondrocytes.

Regents

Tert-Butyl hydroperoxide solution (TBHP) was acquired from Sigma-Aldrich (St Louis, MO, USA). BNIP3 antibody (USA, #44060) was acquired from Cell Signaling Technology. Proteintech provides GAPDH antibody (USA, 60004-1-Ig, diluted 1:1000).

RT-qPCR validation of the hub genes

Total RNA extraction with a total RNA extraction kit (Toyobo, Japan). The First Strand cDNA Synthesis kit (Toyobo, Japan) was used to synthesize complementary DNA (cDNA) from total RNA. Next, the cDNA was amplified with SYBR Green real-time PCR Master Mix (Toyobo, Japan) using the following cycling conditions: 30 seconds of polymerase activation at 95°C followed by 40 cycles of 95°C for 5 seconds and 60°C for 30 seconds. Internal control was glyceraldehyde-3-phosphate dehydrogenase (GAPDH). CDNA samples were run in triplicate. The primer sequences are used in the key resources table. Normalization of consequences to housekeeping gene glyceraldehyde-3-phosphate dehydrogenase (GAPDH) expression using the 2-ΔΔCt algorithm.

Analyses of western blot

Chondrocytes were seeded at a density of 5x10ˆ5 cells/well in six-well plates and allowed to adhere for 24 hours before being subjected to different treatments. The control group remained untreated, while in the experimental group, chondrocytes were exposed to 50 μM TBHP for 24 hours. Western blot analysis was conducted as previously described.54 The cells were washed three times with PBS and lysed using RIPA buffer (Boster, China, AR0102) containing a 1% protease inhibitor cocktail, on ice for 30 minutes. The lysate was centrifuged at 13,000 rpm for 20 minutes to collect the supernatant, and protein concentration was determined using a BCA assay kit (Boster, China, AR0146). A total of 25 μg of protein per sample was subjected to electrophoresis on 12% SDS-PAGE gels and transferred onto PVDF membranes (Millipore, USA). The membranes were blocked with 5% skim milk for one hour at room temperature and incubated with primary antibody overnight at 4°C, followed by secondary antibody for one hour at room temperature. Protein bands were visualized using the ChemiDocTM XRS+ system (Bio-Rad Laboratories, CA, United States) with a chemiluminescent substrate kit (Boster, China), and band density was quantified using ImageJ version 1.48.

Immunofluorescence

A 12-well plate was seeded with chondrocytes. After treatment, the cells were as shown in 4% paraformaldehyde for 15 minutes at room temperature. Following permeabilization with 0.5% TritonX-100 for 20 minutes, the cells were blocked with 5% BSA for 1 hour at room temperature. BNIP3 antibody was then incubated overnight with chondrocytes at 4°C. Cells were treated with goat anti-rabbit secondary antibody for 1 hour in the dark the following day. We incubated the cells with DAPI for 5 minutes after washing them three times with PBS. An Evosflauto fluorescence microscope (Life Technologies, USA) was utilized to obtain the images.

Logistic regression model

Logistic regression is a generalized linear regression analysis model that allows the automatic diagnosis of diseases.55 OA samples were divided into the response variable 1 in this study, while NC samples were represented by the response variable 0. In the first step, we use stepwise regression analysis to eliminate factors that are insignificant for the response variable and thus streamline the model. The Akaike Information Criterion (AIC) is minimized by iteratively adding or removing variables from the model.56,57 These significant factors were then fitted to the response variables for logistic regression. Utilize the receiver operating characteristic curve (ROC) and the area under the ROC curve to evaluate the diagnostic efficacy of the constructed model. "stats" and "pROC" packages were used to perform these analyses.

Immune infiltration and immune-related factors analysis

We employed the CIBERSORT system58 to evaluate immune cell infiltration in the OA microenvironment. This advanced tool, built on linear support vector regression principles, accurately distinguishes among 22 immune cell types, including B cells, T cells, and myeloid cells, by analyzing their specific biomarkers. Our study used gene expression data from the GSE55235 dataset to quantify these cells' relative proportions. We then employed GraphPad Prism to explore the correlation between variations in these immune cells and the identified hub gene, providing insights into immune involvement in OA.

Analysis of protein subcellular localization

The compartment database predicted the subcellular localization of target proteins. Drug-gene interactions were explored using the Drugbank Database and Cytoscape for visualization, while molecular docking involved PubChem and PDB databases for structural analysis.

The compartment database predicted the subcellular localization of target proteins. An overview of the protein's uniform localization can be achieved with a schematic cell.59

Analysis of drug-gene interactions and molecular docking

Drug-gene interactions were explored using the Drugbank Database and Cytoscape for visualization. Drugbank Database60 is used to determine existing and potentially relevant compounds to investigate drug-gene interactions. DrugBank is a repository of chemoinformatics and bioinformatics information on drugs and their targets. Cytoscape was used for data visualization, and molecular structures of ligands and putative proteins were obtained from PubChem,61 and PDB databases.62 Molecular docking involved PubChem and PDB databases for structural analysis. Docking simulations were performed by auto dock vina to generate docking structures, and PyMOL2 software was used for visualization.

Quantification and statistical analysis

We utilized R Software 4.1.1. Pearson correlation analysis was applied to establish relationships between genes. Group comparisons were conducted using the Chi-square test or Fishe R exact test for categorical variables. Continuous variables were analyzed using the independent Student's t-test for normally distributed data, and the Mann-Whitney U test for non-normally distributed datasets. All tests were two-tailed with a significance level set at P < 0.05, ensuring the reliability of our findings. The comprehensive statistical analysis was augmented by GraphPad Prism 7, allowing for a nuanced interpretation of the data. All data represent mean ± S.D. ∗P < 0.05, ∗∗P < 0.01.

Supplemental information

Document S1. Figure S1 and Tables S4 and S6

Table S1. A total of 2,136 differentially expressed genes were identified from GSE55235, related to Figure 2A

This table corresponds to the “Microarrays Data Source and Differential Expression Analysis” section of the STAR Methods, and the differentially expressed genes were obtained under the criteria of |log2FC| > 1.5 and p < 0.05.

Table S2. 269 MElightcyan1 module hub genes, related to Figure 3

The table corresponds to the “Mitochondrial Autophagy-Related Genes and WGCNA” section of the STAR Methods, wherein genes most closely associated with mitochondrial autophagy and OA were identified through WGCNA.

Table S3. 2,364 mitochondrial autophagy-related genes, related to Figure 4

The table pertains to the “Hub Gene Identification and Validation” section of the STAR Methods, and the mitochondrial autophagy-related genes were obtained through the Genecards database.

Table S5. CytoHubba analysis, related to Figure 4

The table is related to the “Hub Gene Identification and Validation” section of the STAR Methods, wherein 31 central genes were filtered down to 6 interacting central genes using cytoHubba.

Supplemental information can be found online at https://doi.org/10.1016/j.isci.2024.110741.
==== Refs
References

1 Park D.R. Kim J. Kim G.M. Lee H. Kim M. Hwang D. Lee H. Kim H.S. Kim W. Park M.C. Osteoclast-associated receptor blockade prevents articular cartilage destruction via chondrocyte apoptosis regulation Nat. Commun. 11 2020 4343 10.1038/s41467-020-18208-y 32859940
2 Story Jovanova O. Nedeljkovic I. Spieler D. Walker R.M. Liu C. Luciano M. Bressler J. Brody J. Drake A.J. Evans K.L. DNA Methylation Signatures of Depressive Symptoms in Middle-aged and Elderly Persons: Meta-analysis of Multiethnic Epigenome-wide Studies JAMA Psychiatr. 75 2018 949 959 10.1001/jamapsychiatry.2018.1725
3 Huang J. Zhao L. Fan Y. Liao L. Ma P.X. Xiao G. Chen D. The microRNAs miR-204 and miR-211 maintain joint homeostasis and protect against osteoarthritis progression Nat. Commun. 10 2019 2876 10.1038/s41467-019-10753-5 31253842
4 Ma Z. Tao C. Sun L. Qi S. Le Y. Wang J. Li C. Liu X. Zhang J. Zhao J. In Situ Forming Injectable Hydrogel For Encapsulation Of Nanoiguratimod And Sustained Release Of Therapeutics Int. J. Nanomed. 14 2019 8725 8738 10.2147/IJN.S214507
5 Yao Q. Wu X. Tao C. Gong W. Chen M. Qu M. Zhong Y. He T. Chen S. Xiao G. Osteoarthritis: pathogenic signaling pathways and therapeutic targets Signal Transduct. Targeted Ther. 8 2023 56 10.1038/s41392-023-01330-w
6 Duong V. Oo W.M. Ding C. Culvenor A.G. Hunter D.J. Evaluation and Treatment of Knee Pain: A Review JAMA 330 2023 1568 1580 10.1001/jama.2023.19675 37874571
7 Duan R. Xie H. Liu Z.Z. The Role of Autophagy in Osteoarthritis Front. Cell Dev. Biol. 8 2020 608388 10.3389/fcell.2020.608388
8 Liu Z. Wang T. Sun X. Nie M. Autophagy and apoptosis: regulatory factors of chondrocyte phenotype transition in osteoarthritis Hum. Cell 36 2023 1326 1335 10.1007/s13577-023-00926-2 37277675
9 Chen P. Zheng L. Wang Y. Tao M. Xie Z. Xia C. Gu C. Chen J. Qiu P. Mei S. Desktop-stereolithography 3D printing of a radially oriented extracellular matrix/mesenchymal stem cell exosome bioink for osteochondral defect regeneration Theranostics 9 2019 2439 2459 10.7150/thno.31017 31131046
10 Alva V. Nam S.Z. Söding J. Lupas A.N. The MPI bioinformatics Toolkit as an integrative platform for advanced protein sequence and structure analysis Nucleic Acids Res. 44 2016 W410 W415 10.1093/nar/gkw348 27131380
11 Liu X. Xiao H. Peng X. Chai Y. Wang S. Wen G. Identification and comprehensive analysis of circRNA-miRNA-mRNA regulatory networks in osteoarthritis Front. Immunol. 13 2022 1050743 10.3389/fimmu.2022.1050743
12 Xia L. Gong N. Identification and verification of ferroptosis-related genes in the synovial tissue of osteoarthritis using bioinformatics analysis Front. Mol. Biosci. 9 2022 992044 10.3389/fmolb.2022.992044
13 Liu H. Deng Z. Yu B. Liu H. Yang Z. Zeng A. Fu M. Identification of SLC3A2 as a Potential Therapeutic Target of Osteoarthritis Involved in Ferroptosis by Integrating Bioinformatics, Clinical Factors and Experiments Cells 11 2022 3430 10.3390/cells11213430 36359826
14 Zhang Y. Zhu T. He F. Chen A.C. Yang H. Zhu X. Identification of Key Genes and Pathways in Osteoarthritis via Bioinformatic Tools: An Updated Analysis Cartilage 13 2021 1457S 1464S 10.1177/19476035211008975 33855867
15 Xu W. Wang X. Liu D. Lin X. Wang B. Xi C. Kong P. Yan J. Identification and validation of hub genes and potential drugs involved in osteoarthritis through bioinformatics analysis Front. Genet. 14 2023 1117713 10.3389/fgene.2023.1117713
16 Zeng J. Jiang X. Jiang M. Cao Y. Jiang Y. Bioinformatics analysis of hub genes as osteoarthritis prognostic biomarkers Sci. Rep. 13 2023 22894 10.1038/s41598-023-48446-1
17 Lin J. Guo Z. Zheng Z. Hou L. Xu J. Liu Q. Du T. Guo F. Jing X. Desferoxamine protects against hemophilic arthropathy through the upregulation of HIF-1α-BNIP3 mediated mitophagy Life Sci. 312 2023 121172 10.1016/j.lfs.2022.121172
18 Zhou Y. Zhou B. Pache L. Chang M. Khodabakhshi A.H. Tanaseichuk O. Benner C. Chanda S.K. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets Nat. Commun. 10 2019 1523 10.1038/s41467-019-09234-6 30944313
19 Hu S. Zhang C. Ni L. Huang C. Chen D. Shi K. Jin H. Zhang K. Li Y. Xie L. Stabilization of HIF-1alpha alleviates osteoarthritis via enhancing mitophagy Cell Death Dis. 11 2020 481 10.1038/s41419-020-2680-0 32587244
20 Zhuo Q. Yang W. Chen J. Wang Y. Metabolic syndrome meets osteoarthritis Nat. Rev. Rheumatol. 8 2012 729 737 10.1038/nrrheum.2012.135 22907293
21 Varela-Eirin M. Loureiro J. Fonseca E. Corrochano S. Caeiro J.R. Collado M. Mayan M.D. Cartilage regeneration and ageing: Targeting cellular plasticity in osteoarthritis Ageing Res. Rev. 42 2018 56 71 10.1016/j.arr.2017.12.006 29258883
22 Kraft C. Peter M. Hofmann K. Selective autophagy: ubiquitin-mediated recognition and beyond Nat. Cell Biol. 12 2010 836 841 10.1038/ncb0910-836 20811356
23 Ansari M.Y. Khan N.M. Ahmad I. Haqqi T.M. Parkin clearance of dysfunctional mitochondria regulates ROS levels and increases survival of human chondrocytes Osteoarthritis Cartilage 26 2018 1087 1097 10.1016/j.joca.2017.07.020 28801211
24 Nakahira K. Haspel J.A. Rathinam V.A.K. Lee S.J. Dolinay T. Lam H.C. Englert J.A. Rabinovitch M. Cernadas M. Kim H.P. Autophagy proteins regulate innate immune responses by inhibiting the release of mitochondrial DNA mediated by the NALP3 inflammasome Nat. Immunol. 12 2011 222 230 10.1038/ni.1980 21151103
25 Tang Q. Zheng G. Feng Z. Chen Y. Lou Y. Wang C. Zhang X. Zhang Y. Xu H. Shang P. Liu H. Trehalose ameliorates oxidative stress-mediated mitochondrial dysfunction and ER stress via selective autophagy stimulation and autophagic flux restoration in osteoarthritis development Cell Death Dis. 8 2017 e3081 10.1038/cddis.2017.453
26 Alvarez-Garcia O. Matsuzaki T. Olmer M. Plate L. Kelly J.W. Lotz M.K. Regulated in Development and DNA Damage Response 1 Deficiency Impairs Autophagy and Mitochondrial Biogenesis in Articular Cartilage and Increases the Severity of Experimental Osteoarthritis Arthritis Rheumatol. 69 2017 1418 1428 10.1002/art.40104 28334504
27 Tian X. Teng J. Chen J. New insights regarding SNARE proteins in autophagosome-lysosome fusion Autophagy 17 2021 2680 2688 10.1080/15548627.2020.1823124 32924745
28 Tahrir F.G. Knezevic T. Gupta M.K. Gordon J. Cheung J.Y. Feldman A.M. Khalili K. Evidence for the Role of BAG3 in Mitochondrial Quality Control in Cardiomyocytes J. Cell. Physiol. 232 2017 797 805 10.1002/jcp.25476 27381181
29 Cen X. Chen Y. Xu X. Wu R. He F. Zhao Q. Sun Q. Yi C. Wu J. Najafov A. Xia H. Pharmacological targeting of MCL-1 promotes mitophagy and improves disease pathologies in an Alzheimer's disease mouse model Nat. Commun. 11 2020 5731 10.1038/s41467-020-19547-6 33184293
30 Lin Q. Li S. Jiang N. Jin H. Shao X. Zhu X. Wu J. Zhang M. Zhang Z. Shen J. Inhibiting NLRP3 inflammasome attenuates apoptosis in contrast-induced acute kidney injury through the upregulation of HIF1A and BNIP3-mediated mitophagy Autophagy 17 2021 2975 2990 10.1080/15548627.2020.1848971 33345685
31 Tu B. Fang R. Zhu Z. Chen G. Peng C. Ning R. Comprehensive analysis of arachidonic acid metabolism-related genes in diagnosis and synovial immune in osteoarthritis: based on bulk and single-cell RNA sequencing data Inflamm. Res. 72 2023 955 970 10.1007/s00011-023-01720-4 36995411
32 Penatti A. Facciotti F. De Matteis R. Larghi P. Paroni M. Murgo A. De Lucia O. Pagani M. Pierannunzii L. Truzzi M. Differences in serum and synovial CD4+ T cells and cytokine profiles to stratify patients with inflammatory osteoarthritis and rheumatoid arthritis Arthritis Res. Ther. 19 2017 103 10.1186/s13075-017-1305-1 28526072
33 Paulissen S.M.J. van Hamburg J.P. Dankers W. Lubberts E. The role and modulation of CCR6+ Th17 cell populations in rheumatoid arthritis Cytokine 74 2015 43 53 10.1038/nri2094 25828206
34 McInnes I.B. Schett G. Cytokines in the pathogenesis of rheumatoid arthritis Nat. Rev. Immunol. 7 2007 429 442 10.3390/biomedicines10092111 17525752
35 Nees T.A. Zhang J.A. Platzer H. Walker T. Reiner T. Tripel E. Moradi B. Rosshirt N. Infiltration Profile of Regulatory T Cells in Osteoarthritis-Related Pain and Disability Biomedicines 10 2022 2111 10.3390/biomedicines10092111 36140212
36 Buckley M.G. Walters C. Wong W.M. Cawley M.I. Ren S. Schwartz L.B. Walls A.F. Mast cell activation in arthritis: detection of alpha- and beta-tryptase, histamine and eosinophil cationic protein in synovial fluid Clin. Sci. 93 1997 363 370 10.1042/cs0930363
37 Park Y.B. Kim J.H. Ha C.W. Lee D.H. Clinical Efficacy of Platelet-Rich Plasma Injection and Its Association With Growth Factors in the Treatment of Mild to Moderate Knee Osteoarthritis: A Randomized Double-Blind Controlled Clinical Trial As Compared With Hyaluronic Acid Am. J. Sports Med. 49 2021 487 496 10.1177/0363546520986867 33523756
38 Faust H.J. Zhang H. Han J. Wolf M.T. Jeon O.H. Sadtler K. Peña A.N. Chung L. Maestas D.R. Jr. Tam A.J. IL-17 and immunologically induced senescence regulate response to injury in osteoarthritis J. Clin. Invest. 130 2020 5493 5507 10.1172/JCI134091 32955487
39 Wang Q. Lepus C.M. Raghu H. Reber L.L. Tsai M.M. Wong H.H. von Kaeppler E. Lingampalli N. Bloom M.S. Hu N. IgE-mediated mast cell activation promotes inflammation and cartilage destruction in osteoarthritis Elife 8 2019 e39905 10.7554/eLife.39905
40 Vane J.R. Botting R.M. The mechanism of action of aspirin Thromb. Res. 110 2003 255 258 10.1016/s0049-3848(03)00379-7 14592543
41 Pietrocola F. Castoldi F. Maiuri M.C. Kroemer G. Aspirin-another caloric-restriction mimetic Autophagy 14 2018 1162 1163 10.1080/15548627.2018.1454810 29929449
42 Wluka A.E. Ding C. Wang Y. Jones G. Urquhart D.M. Cicuttini F.M. Aspirin is associated with reduced cartilage loss in knee osteoarthritis: Data from a cohort study Maturitas 81 2015 394 397 10.1016/j.maturitas.2015.04.015 26027908
43 Yoon J.B. Kim S.J. Hwang S.G. Chang S. Kang S.S. Chun J.S. Non-steroidal anti-inflammatory drugs inhibit nitric oxide-induced apoptosis and dedifferentiation of articular chondrocytes independent of cyclooxygenase activity J. Biol. Chem. 278 2003 15319 15325 10.1074/jbc.M212520200 12588866
44 Munzel T. Gori T. Keaney J.F. Jr. Maack C. Daiber A. Pathophysiological role of oxidative stress in systolic and diastolic heart failure and its therapeutic implications Eur. Heart J. 36 2015 2555 2564 10.1093/eurheartj/ehv305 26142467
45 Xu J. Huang Z. Lin L. Fu M. Song Y. Shen Y. Ren D. Gao Y. Su Y. Zou Y. miRNA-130b is required for the ERK/FOXM1 pathway activation-mediated protective effects of isosorbide dinitrate against mesenchymal stem cell senescence induced by high glucose Int. J. Mol. Med. 35 2015 59 71 10.3892/ijmm.2014.1985 25355277
46 Chang J.W. Liu S.C. Lin Y.Y. He X.Y. Wu Y.S. Su C.M. Tsai C.H. Chen H.T. Fong Y.C. Hu S.L. Nesfatin-1 Stimulates CCL2-dependent Monocyte Migration And M1 Macrophage Polarization: Implications For Rheumatoid Arthritis Therapy Int. J. Biol. Sci. 19 2023 281 293 10.7150/ijbs.77987 36594093
47 Hu X. Ni S. Zhao K. Qian J. Duan Y. Bioinformatics-Led Discovery of Osteoarthritis Biomarkers and Inflammatory Infiltrates Front. Immunol. 13 2022 871008 10.3389/fimmu.2022.871008
48 Cheng Q. Chen X. Wu H. Du Y. Three hematologic/immune system-specific expressed genes are considered as the potential biomarkers for the diagnosis of early rheumatoid arthritis through bioinformatics analysis J. Transl. Med. 19 2021 18 10.1186/s12967-020-02689-y 33407587
49 Qin J. Zhang J. Wu J.J. Ru X. Zhong Q.L. Zhao J.M. Lan N.H. Identification of autophagy-related genes in osteoarthritis articular cartilage and their roles in immune infiltration Front. Immunol. 14 2023 1263988 10.3389/fimmu.2023.1263988
50 Stelzer G. Rosen N. Plaschkes I. Zimmerman S. Twik M. Fishilevich S. Stein T.I. Nudel R. Lieder I. Mazor Y. The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses Curr. Protoc. Bioinformatics 54 2016 1.30.1 1.30.33 10.1002/cpbi.5
51 Meng J. Du H. Lv H. Lu J. Li J. Yao J. Identification of the osteoarthritis signature gene PDK1 by machine learning and its regulatory mechanisms on chondrocyte autophagy and apoptosis Front. Immunol. 13 2022 1072526 10.3389/fimmu.2022.1072526
52 Liberzon A. Birger C. Thorvaldsdóttir H. Ghandi M. Mesirov J.P. Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection Cell Syst. 1 2015 417 425 10.1016/j.cels.2015.12.004 26771021
53 Yao X. Sun K. Yu S. Luo J. Guo J. Lin J. Wang G. Guo Z. Ye Y. Guo F. Chondrocyte ferroptosis contribute to the progression of osteoarthritis J. Orthop. Translat. 27 2021 33 43 10.1016/j.jot.2020.09.006 33376672
54 Yao X. Jing X. Guo J. Sun K. Deng Y. Zhang Y. Guo F. Ye Y. Icariin Protects Bone Marrow Mesenchymal Stem Cells Against Iron Overload Induced Dysfunction Through Mitochondrial Fusion and Fission, PI3K/AKT/mTOR and MAPK Pathways Front. Pharmacol. 10 2019 163 10.3389/fphar.2019.00163 30873034
55 Shen W. Song Z. Zhong X. Huang M. Shen D. Gao P. Qian X. Wang M. He X. Wang T. Sangerbox: A comprehensive, interaction-friendly clinical bioinformatics analysis platform iMeta 1 2022 e36 10.1002/imt2.36 38868713
56 Coat J. Demoersman J. Beuzit S. Cornec D. Devauchelle-Pensec V. Saraux A. Pers J.O. Anti-B lymphocyte immunotherapy is associated with improvement of periodontal status in subjects with rheumatoid arthritis J. Clin. Periodontol. 42 2015 817 823 10.1111/jcpe.12433 26205081
57 Lai D. Tan L. Zuo X. Liu D. Jiao D. Wan G. Lu C. Shen D. Gu X. Prognostic Ferroptosis-Related lncRNA Signatures Associated With Immunotherapy and Chemotherapy Responses in Patients With Stomach Cancer Front. Genet. 12 2021 798612 10.3389/fgene.2021.798612
58 Newman A.M. Liu C.L. Green M.R. Gentles A.J. Feng W. Xu Y. Hoang C.D. Diehn M. Alizadeh A.A. Robust enumeration of cell subsets from tissue expression profiles Nat. Methods 12 2015 453 457 10.1038/nmeth.3337 25822800
59 Binder J.X. Pletscher-Frankild S. Tsafou K. Stolte C. O'Donoghue S.I. Schneider R. Jensen L.J. COMPARTMENTS: unification and visualization of protein subcellular localization evidence Database 2014 2014 bau012 10.1093/database/bau012
60 Wishart D.S. Feunang Y.D. Guo A.C. Lo E.J. Marcu A. Grant J.R. Sajed T. Johnson D. Li C. Sayeeda Z. DrugBank 5.0: a major update to the DrugBank database for 2018 Nucleic Acids Res. 46 2018 D1074 D1082 10.1093/nar/gkx1037 29126136
61 Kim S. Chen J. Cheng T. Gindulyte A. He J. He S. Li Q. Shoemaker B.A. Thiessen P.A. Yu B. PubChem in 2021: new data content and improved web interfaces Nucleic Acids Res. 49 2021 10.1093/nar/gkaa971 D1388-d1395
62 Berman H.M. Westbrook J. Feng Z. Gilliland G. Bhat T.N. Weissig H. Shindyalov I.N. Bourne P.E. The Protein Data Bank Nucleic Acids Res. 28 2000 235 242 10.1093/nar/28.1.235 10592235
