Bioinformatics Analysis and Experimental Validation of Mitochondrial Autophagy Genes in Knee Osteoarthritis

Background

Knee Osteoarthritis (KOA) is a chronic and progressive disease characterized by loss of articular cartilage, subchondral sclerosis, and synovial and periarticular structural abnormalities that cause pain, swelling and even joint deformity in the surrounding tissues of the knee.¹ Its pathogenic factors include abnormal biologic factors, genetic factors, cell aging and apoptosis, local inflammatory factors, free radicals and proteases, etc. Each of these risk factors may cause KOA, KOA is a degenerative, attritional disease,² its pathology mainly involves articular cartilage and synovium, the key factor of KOA cartilage degeneration is chondrocyte death,³ which is caused by apoptosis and autophagy.^4,5 Under normal conditions, articular chondrocytes maintain a dynamic balance between extracellular matrix (ECM) synthesis and degradation,⁶ which is broken in response to injury or inflammation and other stimuli, triggering chondrocyte stress and ECM degradation, it leads to abnormal accumulation of damaged proteins and dysfunction of organelles such as endoplasmic reticulum and mitochondria,^7,8 thus, cartilage degradation, chondrocyte apoptosis, subchondral bone dysfunction, and gradually develop into knee osteoarthritis.⁹ According to the latest Global Burden of Disease study, about 250 million people worldwide suffer from osteoarthritis, and knee joint is the most common site of osteoarthritis.² Due to the lack of early diagnostic indicators, patients with KOA often miss the opportunity for optimal treatment, resulting in a poor prognosis.

Mitochondria are membrane-coated organelles present in eukaryotic cells that regulate important cell functions and cell survival. Mitochondria have many important functions, including cell adenosine triphosphate (ATP) production, Ca2+ buffering, and being the main source of endogenous reactive oxygen species (ROS) under oxidative stress.¹⁰ Mitochondrial autophagy is a major intracellular protective mechanism, which can effectively promote the clearance of damaged mitochondria and thus maintain mitochondrial function.^11,12 Mitochondria can protect chondrocytes and provide energy for cells, which is an important mechanism to maintain cell homeostasis and is closely related to the occurrence and development of cartilage degeneration in osteoarthritis.¹³ Inhibition of autophagy can lead to mitochondrial dysfunction and catabolic imbalance, leading to chondrocyte injury and apoptosis,¹⁴ thus leading to osteoarthritis.

In recent years, bioinformatics has provided powerful strategies for the screening of molecular markers, and cell classification tools have also facilitated the analysis of infiltration patterns of disease immune cells. Based on the role of mitochondrial autophagy in osteoarthritis, we explored the mechanism of mitochondrial autophagy gene in knee arthritis, attempted to identify diagnostic markers of KOA, and further explored the role of immune cell infiltration in KOA, providing theoretical possibilities for the development of new targeted drugs, which is of great significance for the study of the pathogenesis and prevention of KOA.

Materials and Methods

Data Extraction

RNA-sequencing (RNA-seq) data of KOA patients (GSE114007, GSE169077 and GSE51588) were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE114007 dataset, including cartilaginous tissues from 20 KOA samples and 18 normal samples, was used as the training set. The GSE169077 dataset, containing 6 KOA samples and 5 normal samples, was used as the validation set. The GSE51588 dataset (KOA: normal = 40:10) was also used as validation set. In addition, 88 mitochondrial autophagy gene (MGs) were downloaded for this study based on the relevant literature (Sun et al 2021).

Differential Expression Analysis

In the GSE114007 dataset, differentially expressed genes (DEGs) between KOA samples and normal samples were detected utilizing the “DESeq2” R package (version 1.34.0) with an adjusted p<0.05 and a |log2FC|>0.5. Clinical characteristics (gender and age) were also controlled for to exclude the effect of other confounding factors. Subsequently, the DEGs were intersected with MGs using the “VennDiagram” R package (version 1.6.20) (Xin et al 2022) to get differentially expressed MGs (DEMGs). Subsequently, The Gene Ontology (GO) and The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of the DEMGs were performed using the “clusterProfiler” R package (Jahr et al 2019; Lazarou 2015).

PPI Network Construction and Hub DEMGs Screening

To explore the protein interactions among DEMGs, the STRING (https://string-db.org) website was used to construct a PPI network, and DEMGs with interaction scores greater than 0.4 were retained. Moreover, hub DEMGs were identified using MCODE algorithm in Cytoscape (Xin et al 2022).

Construction and Validation of a Diagnostic Model for DEMGs

Based on the identified hub DEMGs, the SVM-RFE method was applied to screened the model DEMGs using the “caret” R package (version 6.0–91, https://CRAN.R-project.org/package=caret). Subsequently, diagnostic value of the model DEMGs was validated in a validation set (GSE169077). In addition, we used the “pROC” R package (version 1.17.0.1) to map the ROC curve to assess their diagnostic value. Finally, the model DEMGs with diagnostic value were validated in the training and validation set.

Immune Cell Infiltration with CIBERSORT

To further explore the correlation between model DEMGs and different immune cells, the infiltration of 22 immune cells in the GSE114007 dataset was estimated by CIBERSORT algorithm. Then, the “psych” R package (version 2.2.3) was performed between model DEMGs and immune cells.

Experimental Verification

Cartilage tissue samples were obtained from KOA patients and normal samples with their knowledge and consent from Beijing Jishuitan Hospital Guizhou Hospital, and this study was approved by the Beijing Jishuitan Hospital Guizhou Hospital ethics committee, all methods were performed in accordance with relevant guidelines and regulations in the《Ethical considerations in the review of biomedical research involving human subjects》and the《Declaration of Helsinki》(approval number: LW20221206).

10 pairs of frozen cartilage tissue samples were divided into two groups, of which 10 samples were normal group (Normal) and the other 10 samples were KOA group (Case). Then, total RNA of samples was isolated and purified by TRIzol (Ambion) reagent following the instruction manual. Then, the extracted RNA was tested for concentration by NanoPhotometer N50. Next, reverse transcription via SureScript-First-strand-cDNA-synthesis-kit (Servicebio) by an ordinary PCR instrument. Reverse transcription product cDNA was diluted 5–20 times with ddH₂O (RNase / DNase free). Subsequently, polymerase chain reaction (PCR) amplification reaction was performed by CFX96 real-time quantitative PCR instrument. 1 min at 95 °C (pre-denaturation), followed by at 95 °C for 20s (denaturation), 55 °C for 20s (annealing) and 72 °C for 30s (elongation). The above reactions were subjected to forty cycles. Primer sequences were showed in Table 1.

Table 1 The Primer Sequences of Diagnostic Genes

Statistical Analysis

All statistical analyses and visual plotting of the results were performed based on R software (https://www.r-project.org/, version 4.0.3, R Statistical Computing Project). Correlation of model DEMGs with immune cells using Spearman coefficient analysis.

Results

Identification of DEMGs

To screen DEMGs in KOA, we firstly screened DEGs between KOA samples and normal samples in the GSE114007 dataset. A total of 4279 DEGs, including 1969 down-regulated and 2310 up-regulated in KOA samples, were identified (Figure 1a and b). Comparison of differences in clinical characteristics revealed no significant effect. (Table S1). Then, 15 DEMGs were obtained for subsequent analysis by intersecting DEGs and MGs (Figure 1c).

Figure 1 Differential expression analysis. (a) Differentially expressed gene expression volcano plot between KOA and normal controls. (b) Differentially expressed gene expression heatmap between KOA and normal controls (adjusted p<0.05 and a |log2FC|>0.5). (c) Venn diagram shows the intersection of DEGs and MGs. Abbreviations: KOA, knee osteoarthritis; FC, fold change; DEGs, differentially expressed genes; MGs, mitochondrial autophagy gene.

GO and KEGG Analysis of the DEMGs

To explore the regulatory pathways and biological functions associated with the expression of DEMGs, we performed GO and KEGG analyses on these DMAEGs in the GSE114007 dataset. GO enrichment analysis suggested that in biological process (BP), these DEMGs were mainly involved in 5 items such as organelle disassembly, macroautophagy and cellular component disassembly; these DEMGs were mainly involved in 5 items such as mitochondrial outer membrane, organelle outer membrane and outer membrane in cellular component (CC); in molecular function (MF), these DEMGs are mainly involved in two items, ubiquitin protein ligase binding and ubiquitin-like protein ligase binding (Figure 2a). Moreover, KEGG pathway enrichment analysis indicated that these DEMGs were significantly associated with Mitophagy-animal, Shigellosis, Autophagy-animal, and FoxO signaling pathway pathways (Figure 2b).

Figure 2 Functional enrichment analysis. (a) GO enrichment analysis based on differentially expressed genes was performed in the GSE114007 databases. (b) KEGG enrichment analysis was performed to explore the associated pathways in the GSE114007 databases. Abbreviations: GO, Gene Ontology; KEGG, The Kyoto Encyclopedia of Genes and Genomes; BP, biological progress; CC, cellular component; MF, molecular function.

Identification of the Hub DEMGs Through PPI Network

To explore the interaction between 15 DEMGs, the STRING website was applied to build a PPI network. Ultimately, a PPI network including 36 interactions and 13 nodes was obtained (Figure 3a). In addition, 8 hub DEMGs (ULK1, CALCOCO2, MAP1LC3B, BNIP3L, GABARAPL1, BNIP3, FKBP8, and FOXO3) were screened by the MCODE algorithm (Figure 3b).

Figure 3 The PPI network of DEMGs. (a) The PPI network shows the DEMGs with a minimum required interaction score >0.4. (b) The most significant module was obtained from PPI network with 8 hub genes. Abbreviations: PPI, protein-protein interaction; DEMGs, differentially expressed mitochondrial autophagy gene.

Construction and Validation the DEMGs Diagnostic Model for KOA Patients

To further validation and selection of the DEMGs with significantly characteristic value of classifying KOA and normal samples from GSE114007 dataset, the SVM-RFE algorithm was performed to identify 5 model DEMGs (BNIP3L, BNIP3, MAP1LC3B, ULK1, and FOXO3) (Figure 4a). Meanwhile, we evaluated and validated the diagnostic value of the 5 model DEMGs in the GSE114007 and GSE169077 datasets, respectively. As shown in Figure 4b, the AUC values of these genes were above 0.9 in the training set, indicating that these genes have a strong diagnostic value. However, in the validation set, only BNIP3 and FOXO3 had AUC values, and its AUC values greater than 0.7 (Figure 4c). Additionally significant differences and consistent trends in BNIP3 and FOXO3 expression were found in the training and validation sets by expression analysis (Figure S1). Thus, these two genes were identified as having strong diagnostic value and were used for subsequent analysis.

Figure 4 Identification of model DEMGs. (a) The SVM-RFE algorithm was performed to identify 5 model DEMGs (BNIP3L, BNIP3, MAP1LC3B, ULK1, and FOXO3). (b) ROC analysis of 5 model DEMGs on the training set. (c) ROC analysis of 5 model DEMGs on the validation set. Abbreviations: SVM-RFE, support vector machine recursive feature elimination; ROC, receiver operating characteristic; AUC, area under the curve.

Correlation Analysis of BNIP3 and FOXO3 in Immune Cells

To accurately assess the composition of immune cells in the KOA tissue microenvironment, the CIBERSORT algorithm was used to estimate the abundance of immune infiltrating cells of 22 different cell types and functional states between KOA and normal samples in the GSE114007 dataset. The results revealed that of the 22 immune cells only naive B cells, M0 Macrophages, and M2 Macrophages were significantly different between the KOA and normal samples (Figure 5a). Furthermore, the correlation analysis of these immune cells with BNIP3 as well as FOXO3 showed that BNIP3 was significantly negatively correlated with M0 Macrophages (|cor|>0.3 and p value<0.05); while FOXO3 was significantly negatively correlated with M0 Macrophages and showed a significant positive correlation with naive B cells (Figure 5b and c).

Figure 5 Immune infiltration analysis. (a) Violin diagram of the proportion of 22 different immune cell types and functional states between KOA and normal samples. (b) Correlation between BNIP3 and immune cells. (c) Correlation between FOXO3 and immune cells. *p < 0.05. Abbreviation: NS, not significant.

Quantitative Real-Time PCR (qPCR) Identification

Based on the qPCR verification results, it can be seen that FOXO3 and BNIP3 were all down-regulated in KOA (p < 0.01), and the validation results are consistent with the above analysis (Figure 6, Table 2).

Table 2 qPCR Verification Results

Figure 6 qPCR validation of the BNIP3 and FOXO3 between KOA (Case) and normal controls. **p < 0.01, ***p < 0.001. Abbreviation: NS, not significant.

Discussion

Mitochondrial autophagy can protect human chondrocytes from a series of pathological effects caused by mitochondrial dysfunction and effectively slow down the progression of cartilage degeneration in KOA. In the past, it was thought that articular chondrocytes were in an ischemic and anoxic environment, and the role of mitochondria in the metabolic activities of chondrocytes was relatively limited. However, recent studies have shown that chondrocytes on the surface of cartilage are in a relatively aerobic state, and mitochondria-mediated aerobic respiration is still an important energy supply mode for articular chondrocytes.¹⁵ Mitochondria are involved in catabolism and anabolism, production of ROS, cell apoptosis and signal transduction. When the body is subjected to oxidative stress, endoplasmic reticulum stress, hypoxia ischemia, nutrient deficiency and other injuries, Mitochondrial autophagy disorders lead to mitochondrial dysfunction and catabolic imbalance.¹⁶ Mitochondrial damage can lead to the release of interleukin-1β (IL-1β), which further leads to chondrocyte death and ECM degradation. Overactivation of mitochondrial autophagy can lead to autophagy death and osteoarthritis.^14,17

In this study, a total of 4279 DEGs were obtained from the GSE114007 dataset, 15 DEMGs were obtained by crossing DEGs and MGs. GO enrichment analysis showed that these DEMGs were mainly involved in organelle decomposition, macroautophagy and cell component decomposition in biological processes (BP). In the same type of cell, autophagy can be protective or harmful, depending on the type and degree of stimulation. Maintaining a healthy mitochondrial population is essential for cell survival. As a kind of autophagy, macrophage plays an important and complex role in bone metabolism. Eukaryotic cells can regulate molecular degradation and organelle renewal through macrophage, and play an important role in bone regeneration, bone metabolism, and bone degenerative diseases such as osteoarthritis.¹⁸ In addition, the main BP of KOA included inflammatory response regulation, collagen catabolic process, extracellular matrix decomposition, etc. When mechanical injury or inflammatory stimulation exceeds the compensatory capacity of articular cartilage, chondrocyte stress and ECM degradation are triggered, thus causing cartilage degeneration and gradually developing into KOA.^19,20 KOA can be alleviated by inhibiting inflammation, chondro breakdown, and extracellular matrix degradation.²¹ In the cell component (CC), these DEMGs mainly involved 5 items such as mitochondrial outer membrane, organelle outer membrane and outer membrane. Mitochondrial outer membrane is a biological membrane located in the outermost layer of mitochondria. Under stimulation, Mitochondrial outer membrane permeabilization (MOMP) can permeabilization. ROS released by damaged chondrocytes promote the expression of matrix metalloproteinases such as MMP13. Meanwhile, excessive ROS generation leads to mitochondrial dysfunction and enhanced permeability of mitochondrial outer membrane.²² Mitochondrial proteins are released into the cytoplasm to mediate downstream reactions, Bnip3 is localized in the mitochondrial outer membrane, it induces mitochondrial autophagy in the presence of mitochondrial membrane permeability and Bax/Bak deficiency.²³ In molecular function (MF), these DEMGs are mainly involved in ubiquitin-protein ligase binding and ubiquitin-like protein ligase binding. The ubiquitin-proteasome system is an important regulatory pathway to maintain protein homeostasis and function in cells, and abnormal regulation of ubiquitination can lead to metabolic bone disease. E3 ubiquitin ligases (E3) are key regulatory factors in the proliferation and differentiation of osteoblasts, play an important role in the regulation of bone formation related proteins and bone turnover, and have been shown to be ideal therapeutic targets for promoting bone formation and reducing bone loss,²⁴ The degradation of BNIP3 by ubiquitination under severe hypoxia condition can lead to the inhibition of mitochondrial autophagy.²⁵ KEGG pathway enrichment analysis showed that these DEMGs are mainly related to mitochondrial and autophagy signaling pathways, FoxO signaling pathway plays a function regulation role in bone cells, and is associated with bone metabolic diseases such as osteoarthritis and osteoporosis.²⁶ FoxO transcription factor can prevent cellular and biological senescence, and its expression in OA cartilage is decreased with aging.²⁷ Promote cartilage regeneration by up-regulating FOXO signaling pathway and alleviate the progression of osteoarthritis.²⁸

In order to explore the interaction between DEMGs, we constructed a PPI network with the differential expression of mitochondrial autophagy genes using STRING, and used the MCODE algorithm to screen hub genes. Then, the SVM-RFE algorithm was used to verify the diagnostic markers of KOA, and the diagnostic value was evaluated and verified. Finally, BNIP3 and FOXO3 have been identified as potential biomarkers for the treatment of KOA patients, play a key role in the progression of KOA and can be used as therapeutic targets. BNIP3 is a mitochondrial outer membrane protein belonging to the Bcl-2 protein family, which regulates apoptosis, autophagy, mitochondrial transformation and cell protection.²⁹ BNIP3 is an important mediator of autophagy.³⁰ As a dual regulator, it has pro-apoptotic function.³¹ and controversial anti-apoptotic function.^32,33 BNIP3 induces mitochondrial autophagy under various stress states. For example, under hypoxia conditions, BNIP3 promotes mitochondrial depolarization and directly activates mitochondrial autophagy by transporting mitochondria to autophagosomes.^34,35 Bnip3 knockout reduces mitochondrial autophagy, resulting in accumulation of damaged mitochondria and increased production of reactive oxygen species.³² BNIP3 induced mitochondrial dysfunction mainly by activating Bax/Bak and the mitochondrial permeability transition pore (mPTP).³⁶ Overexpression of BNIP3 inhibits HIF-1α knockout from enhancing H/R-induced apoptosis and ROS production, reverses the decline in mitochondrial autophagy during injury,³⁷ and protects mitochondrial structure. The loss of BNIP3 leads to the accumulation of abnormal mitochondria,³⁸ aggravates the damage of mitochondrial structure, significantly reduces mitochondrial autophagy, and thus promotes apoptosis.³⁹ FOXO3 is a member of the Fox O transcription factor family, regulating cell survival, apoptosis, autophagy, oxidative stress and inflammatory response.^40,41 It is a key transcription factor involved in pro-inflammatory cytokine induced chondrocyte metabolic disorders, and overexpression prevents IL-1β or TNF-α induced inflammation in various cell types.^42,43 FOXO3 regulates apoptosis of chondrocytes⁴⁴ and maintains meniscus and cartilage homeostasis,^27,45 which may play an important role in the development of KOA. FOXO3 knockout mice resulted in more severe and earlier age-related OA-like changes at 18 months, and the primary cause appeared to be a decrease in protective genes (autophagy, antioxidants).²⁷ Previous research found that FoxO1 can activate mitochondrial autophagy by activating the downstream target gene BNIP3.⁴⁶ FOXO3 increases the expression of BNIP3 by binding to the upstream promoter region of BNIP3, thereby regulating mitochondrial function and integrity.⁴⁷ In order to verify the expression of these autophagy related genes in KOA, we took clinical normal cartilage and KOA cartilage for qPCR verification, and the results showed that FOXO3 and BNIP3 were down-regulated in KOA cartilage, and the verification results were consistent with the above analysis. BNIP3 and FOXO3 were identified as biomarkers for the diagnosis of KOA, which may provide new insights into the pathogenesis of KOA.

The immune system plays a very important role in the development of osteoarthritis. Immune cell infiltration mediates the autoimmune response of osteoarthritis, induces the secretion of chemokines, pro-inflammatory cytokines and proteases, and thus disturbs the immune balance to accelerate cartilage erosion.^48–50 To find specific diagnostic markers and analyze the infiltration pattern of KOA immune cells is of great significance for improving the prognosis of KOA patients. The results showed that naive B cells, M0 Macrophages and M2 Macrophages differ significantly between KOA and normal samples. The proportion of naive B cells, T cells CD8 and follicular helpers T-cells in the articular cartilage tissue of KOA patients was less. The memory resting T-cell CD4, M0 Macrophages and M2 Macrophages accounted for a larger proportion. Correlation analysis between BNIP3 and FOXO3 showed that BNIP3 and FOXO3 were negatively correlated with memory resting T-cell CD4, M0 Macrophages and M2 Macrophages. BNIP3 was positively correlated with T-cell CD8 and FOXO3 was positively correlated with B cells naive. M0 Macrophage was stimulated by interferon IFN-γ and lipopolysaccharide LPS to obtain M1-type macrophages and secrete a large number of pro-inflammatory factors. M0 Macrophage was stimulated by interleukin4 (IL-4) to obtain M2-type macrophages, which had anti-inflammatory and tissue repair effects.⁵¹ We speculate that BNIP3 and FOXO3 affect T-cell CD4 memory resting, M0 Macrophages and M2 Macrophages, and thus participate in the occurrence and progression of KOA. Previous studies have found that the inflammation of KOA is characterized by immune cell infiltration and cytokine secretion, but different from the inflammation of RA, the content of macrophages, T cells and B cells in the synovium of patients is higher than that of healthy patients but lower than that of RA patients.⁵² End-stage KOA articular-derived T cells polarize into CD3 + CD4 + CD8-T cell subsets, increasing the polarization of CD4 + T cells towards activated Th1 cells and cytokine secretion, thereby increasing local inflammation.⁵³ The specific effects of these differentially expressed chemokines on immune infiltration of KOA remain to be further studied.

Conclusion

In this study, we found that BNIP3 and FOXO3 were diagnostic biomarkers of KOA. We also found that naive B cells, M0 Macrophages, and M2 Macrophages may play a key role in the occurrence and development of KOA, and further exploration of these immune cells may determine the target of KOA immunotherapy. However, our study also has certain limitations, most of the current studies remain theoretical, and the relationship between key genes and immune infiltration as well as their role in the development of KOA needs further research. In order to achieve the ultimate goal of KOA treatment, further clinical verification is needed.