Back to Journals » Journal of Inflammation Research » Volume 15

Identification of Key Genes and Potential Mechanisms Based on the Autophagy Regulatory Network in Osteoclasts Using a Murine Osteoarthritis Model

Authors Zhang H, Sun H, Zhang W, Xu Y, Geng D

Received 5 January 2022

Accepted for publication 29 March 2022

Published 12 April 2022 Volume 2022:15 Pages 2333—2347

DOI https://doi.org/10.2147/JIR.S354824

Checked for plagiarism Yes

Review by Single anonymous peer review

Peer reviewer comments 3

Editor who approved publication: Dr Monika Sharma



Video abstract of "Molecular mechanisms of autophagy in osteoclasts involved in OA" [ID 354824].

Views: 573

Haifeng Zhang,* Houyi Sun,* Wei Zhang,* Yaozeng Xu, Dechun Geng

Department of Orthopedics, the First Affiliated Hospital of Soochow University, Suzhou City, People’s Republic of China

*These authors contributed equally to this work

Correspondence: Dechun Geng; Yaozeng Xu, Department of Orthopedics, the First Affiliated Hospital of Soochow University, Suzhou City, People’s Republic of China, Tel +86 512-67780999 ; +86 512-67780249, Email [email protected]; [email protected]

Background: Osteoarthritis (OA) is a degenerative joint disease that acts as a major cause of early disability in the old population. However, the molecular mechanisms of autophagy in osteoclasts involved in OA remain unclear.
Methods: The gene expression profiles were downloaded from the Gene Expression Omnibus (GEO) repository. The NCBI GEO2R and ScanGEO analysis tool were used to identify differentially expressed genes (DEGs). The protein-protein interaction (PPI) network was predicted by the STRING website and visualized with Cytoscape software. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway were performed to enrich GO terms and signaling pathways using Metascape database. To predict LC3-interacting region (LIR) motif among these DEGs, the iLIR database was selected to assess specific short linear sequences. To obtain potential upstream miRNA targets of these DEGs, the mRNA-miRNA interaction networks were predicted by miRWalk database. The knee OA model was performed in mice, and autophagy related mRNAs of osteoclasts were identified. Experimental specimens were further verified with histopathological staining.
Results: Becn1, Atg3, Atg12, Pik3c3, and Gabarapl2 were obtained as coexpressed differential genes. PPI network was constructed and deduced the other 60 related genes. GO and KEGG enrichment networks indicated that autophagy-animal, selective autophagy, and mitophagy mainly participated in autophagy regulation in osteoclasts. The possible LIR sequences were collected to predict motifs. The mRNA–miRNA interaction networks suggested that many miRNAs could regulate autophagy-related genes individually and collectively. The RT–PCR results suggested that these five genes were upregulated in the OA group. Histopathological staining revealed that osteoclasts were increased in subchondral bone, and higher expression of these DEGs in the OA group was compared to the sham group.
Conclusion: Our results reveal that the role of autophagy in osteoclasts could be a regulatory mechanism in OA and that these autophagy-related genes might be targets for the intervention of OA disease.

Keywords: osteoclasts, autophagy, osteoarthritis, bioinformatics, LIR, miRNA

Introduction

Osteoarthritis (OA) is a common degenerative disease and affects the whole joint, including progressive degeneration of articular cartilage, subchondral bone remodeling, osteophyte formation, meniscus degenerative changes, the infrapatellar fat pad and synovial inflammation.1 Worldwide, it is estimated that OA disease implicates 250 million people currently and cause significant health, socioeconomic problems.2 Despite many risk factors have been identified, including aging, genetic disorders, mechanical abnormalities and obesity, the exact etiology of OA remains unclear.3

Osteoclasts are multinucleated giant cells derived from hematopoietic osteoclast precursors and are key regulators of bone resorption.4 Abnormal osteoclast activation is always involved in arthritic bone diseases.5 An increased number and overactivated function of osteoclasts in subchondral bone plates have been observed in OA patients.6 As a result, the thickness of the subchondral bone plate decreases, and the remodeling rate increases, especially in the early stage of OA disease.7 Other results also indicate that abnormally activated osteoclasts could be an important origin of OA pain.8 The activation of osteoclasts play a significant role in initiating and accelerating the pathophysiological process of osteoarthritis.9

Autophagy is a highly conserved catabolic process in eukaryotic cells. It degrades aging or damaged organelles, breaks down unnecessary macromolecules or pathogens, and releases nutrients and energy through lysosomal mechanisms to guarantee cell survival and maintain cell homeostasis.10 Deregulation of autophagy is associated with the pathogenesis of various diseases, including immune disorders, neurodegenerative diseases, cancer and aging.11 Recent studies have shown that autophagy and its related proteins are involved in the process of activating of osteoclasts, such as osteoclastic differentiation, osteoclast-mediated bone resorption and the migration of osteoclasts by regulating different molecular mechanisms.12–14 In inflammatory arthritis, osteoclast precursor cells also activate autophagy to aggravate the development of disease.15

Many autophagy-related genes and signaling pathways were interconnected with the OA process, including AMPK, ULK1 and mTOR.16 These findings would be helpful for understanding the role of autophagy in OA, particularly in cartilage cell research areas. However, the potential autophagy mechanism of osteoclast activation in OA is still vague, and no studies have fulfilled an integrated analysis to identify related biomarkers or gene networks, which represent potential regulatory targets of OA.

To elucidate the autophagy function of osteoclasts involved in OA disease, bioinformatics methods were selected to screen DEGs and predict enrichment networks as well as molecular interactive networks. In addition, the LC3-interacting region sequence and upstream miRNAs of DEGs were identified. And then, autophagy-related differentially expressed genes were verified by experimental OA models in mice. Our study would firstly identify the integrated autophagy regulatory network in osteoclasts and offers new options for exploring the pathogenesis of osteoarthritis.

Methods

Microarray Gene Data Collection and Data Processing

From the GEO Database (https://www.ncbi.nlm.nih.gov/geo/), our team searched and screened the microarray gene datasets with the keywords “osteoarthritis” and “osteoclasts” following the main points: (1) detailed information and reliable source; (2) expression profiling by array; (3) organism came from Mus musculus. According to the above filter direction, GEO accessions GSE41342 and GSE57468 were selected for analysis.

We then compared knee joint tissue including articular cartilage, subchondral bone, meniscus, and the joint capsule with synovium from the sham group versus the DMM model group expression gene dataset (GEO accession GSE41342) and bone marrow-derived macrophage cells (BMMs) induced by RANKL versus BMMs not induced by the RANKL expression dataset (GEO accession GSE57468) in the GEO database. The two microarray datasets were from the following platforms: GPL1261: [Mouse430_2] Affymetrix Mouse Genome 430 2.0 Array; GPL6885: Illumina MouseRef-8 v2.0 expression beadchip.

For the analysis of GSE41342 microarray data, the DEGs were identified using the NCBI GEO2R online tool (https://www.ncbi.nlm.nih.gov/geo/geo2r/). We compared the sham group and OA group at 4 weeks after surgery. Benjamini & Hochberg correction was selected an alternative P-value adjustment method and P < 0.05 was chose as cut-off value. As an online R language analysis tool and has been used to identify pathways of interest, the ScanGEO website (http://scangeo.dartmouth.edu/ScanGEO/) tool was selected to analyze the GSE57468 database.17 The KEGG pathway was chosen using the keyword “autophagy”. To obtain autophagy-related differential genes, we set up a significance level of variance (P < 0.05). Then, the intersection of DEGs between the above two methods was chosen for further analysis of genes.

Construction of the Protein–Protein Interaction Network

The STRING 11.0 (https://string-db.org/) database is a convenient analysis tool to predict and construct a protein–protein interaction network.18 To explore and visualize the biological relationships based on DEGs, we set the parameters to default values and selected a confidence score > 0.4 as the cut-off criterion. Then, Cytoscape software (http://www.cytoscape.org/) was selected to visualize the interaction network among these proteins.19 To finish precise prediction of selected proteins according to weighted degree, the cytoHubba plug-in was used by maximal clique centrality (MCC) algorithm.

Enrichment Analysis and Construction of the Interacted Network

Metascape (http://metascape.org) is gene function enrichment and classification tool used to analyze the main function of the gene sets.20 Gene ontology (GO) included three main components: biological process (BP), cellular component (CC), and molecular function (MF), which is used to annotate gene functions and interconnections. KEGG is an extensive database that maps the pathway annotation results. GO function and KEGG pathway enrichment networks were performed to analyze genes on the Metascape website. P value < 0.05 as cutoff value, Min Enrichment name not less than 3 and Enrichment factor > 3.0 were collected and grouped into clusters. Then we selected the MCODE algorithm to identify interconnected hub genes. Each MCODE network is assigned a unique color representing closely interacting relationships between molecules. In the process of interaction enrichment analysis, min and max network size were selected 3 and 500, respectively.

Identification of LIR Motifs

To distinguish whether these differential proteins are LC3-interacting region containing proteins (LIRCPs) or Atg8-interacting proteins and predict short linear interacting motifs with Atg8-family proteins, the iLIR online database (https://ilir.warwick.ac.uk/) was selected to mine information.21 First of all, we acquired the protein sequence in FASTA format from the NCBI website and then searched the LIR motif in the iLIR autophagy database. From the start to end position, sequences of the LIRCPs as well as the PSSM scores were downloaded and summarized.

Prediction miRNAs and Construction miRNA-mRNA Regulatory Network

The miRWalk 2.0 online web tool (http://mirwalk.umm.uni-heidelberg.de/), as an integrated database that provides predictions of miRNA-target interactions,22 was selected to forecast upstream miRNAs of differential genes and then construct an interactive network of miRNAs–mRNAs. miRNA data targeting the 3ʹUTR region of DEGs were obtained, and the top 20 interacting nodes of each gene were analyzed using CytoHubba plug-in.

Construction of OA Model and Cell Culture

All procedures and experiments were approved by the Animal Ethics Committee of the First Affiliated Hospital of Soochow University and complied with the guidelines of the Care and Use of Laboratory Animals. Eight-week-old male C57BL/6 mice (provided by the Experimental Animal Center of Soochow University) were selected to induce OA model according to previous description.23 The mice were operated for destabilization of the medial meniscus (DMM) in the right knee in the OA group and the contralateral knee joint was regarded as the sham group. Then the mice were euthanized with carbon dioxide at 4 weeks postoperation (n = 6 per group).

Bone marrow-derived macrophage cells (BMMs) as osteoclast precursors were extracted from the marrow of tibia and femur in the sham and OA group. Bone marrow cavity cells were washed and centrifuged, and then cultured with DMEM (HyClone, UT, USA) containing 10% FBS (Gibco, USA). After 16h, nonadherent cells were harvested and incubated with DMEM containing 10% FBS and 30ng/mL M-CSF (R&D Systems, MN, USA) for 3 days. BMMs were the adherent cells in culture plate and used for the follow-up experiment.

Experimental Verification via RT-PCR and Histological Observation

To establish the model in mice, the knee joint samples were harvested, fixed and decalcified at 4 weeks post operation. Then these samples were dehydrated, embedded and sectioned. Hematoxylin and eosin (HE) staining, Safranin O-Fast Green staining along with Osteoarthritis Research Society International (OARSI) scoring were performed to evaluate histopathological difference in two groups conventionally. Tartrate-resistant acid phosphatase (TRAP) staining was applied to assess the activity of osteoclasts. According to the instruction in TRAP Staining Kit (Suzhou Bizhong Biotechnology, Suzhou, China), positive osteoclasts could be stained step by step.

To verify the differential gene expression, total RNA was extracted with standard protocols when BMMs were simulated by 30ng/mL M-CSF and 50ng/mL RANKL (R&D Systems, MN, USA) for 5 days. Real-time polymerase chain reaction (RT–PCR) was used to evaluate the mRNA expression of Becn1, Atg3, Atg12, Pik3c3 and Gabarapl2 in the two groups. The primer sequences were as follows: forward: 5′-TAATAGCTTCACTCTGATCGGG-3′ and reverse: 5′-CAAACAGCGTTTGTAGTTCTGA-3′ for Becn1; forward: 5′-TACGACCTGTACATCACTTACG-3′ and reverse: 5′-GAGGATGGTTTTCAATCAC-3′ for Atg3; forward: 5′-GCCTCGGAACAGTTGTTTATTT-3′ and reverse: 5′-CAGTTTACCATCACTGCCAAAA-3′ for Atg12; forward: 5′-CATGCTCCGACCTCTATGTGACTTG-3′ and reverse: 5′-TTCAGCCACTCGTTCCAATTCCATC-3′ for Pik3c3; forward: 5′-CATCTTCCTGTTTGTGGACAAG-3′ and reverse: 5′-CAGAAGCCAAAAGTGTTCTCTC-3′ for Gabarapl2; forward: 5′-GGTTGTCTCCTGCGACTTCA-3′ and reverse: 5′-TGGTCCAGGGTTTCTTACTCC-3′ for GAPDH; All primers were provided by Sangon Biotech (Shanghai, China).

In order to detect the expression of autophagy-related genes between the sham group and OA group, immunohistochemistry (IHC) staining and immunofluorescence (IF) staining were performed. Briefly, the samples were fixed, decalcified, embedded and sliced. Then, the sections were dewaxed, followed by antigen repair and then blocked with horse serum. Next, these specimens were incubated with primary antibodies and their corresponding secondary antibodies or fluorescent-labeled secondary antibodies. The protein expression of Becn1 (A7353, ABclonal, Wuhan, China), Atg3 (A5809, ABclonal, Wuhan, China), Atg12 (ab155589, Abcam), Pik3c3 (A4021, ABclonal), and Gabarapl2 (A7782, ABclonal) were evaluated though IHC staining. The protein expression and localization of Becn1 were analyzed with IF staining. All sections were observed under an AxioCam HRC microscope (Carl Zeiss, Jena, Germany). Among these, the fluorescence relative intensity of Becn1 was analyzed by Image-Pro Plus 6.0.

Statistical Analysis

The data were processed using GraphPad Prism 8 software (GraphPad Software Inc., La Jolla, CA) and presented as means ± standard deviations. The data were analyzed by Student’s t-tests for comparison between two groups. For the OARSI scores, the Mann–Whitney U-test was used. P < 0.05 was considered significant (*P < 0.05).

Results

Identification of DEGs

To screen the potential DEGs, GSE41342 and GSE57468 microarray data from GEO database were collected and analyzed. With the help of GEO2R and ScanGEO online analysis tools jointly, a total of 34 autophagy-related DEGs were found, among which five genes, such as Becn1, Atg3, Atg12, Pik3c3 and Gabarapl2, were obtained as coexpressed differential genes (Table 1).

Table 1 The Screened DEGs Between Two Datasets

Construction and Analysis of the Protein–Protein Interaction Network

In order to predict interactive candidate genes of these five key genes, a PPI network was built in the STRING database. Five key genes as key nodes were imported, and the other 60 genes were deduced. The PPI enrichment network consisted of 65 nodes and 834 interactions lines (Supplementary Table 1). The average node degree was 25.7 and the average local clustering coefficient was 0.773, indicating it was a satisfactory multicomponent interaction, with a P value <1.0e-16.

To screen important and valuable key genes, the interaction data and node degree of all 65 genes were entered into Cytoscape software. According to the score values of these genes, CytoHubba plug-in could calculate and then rank the gene list by different node sizes and circle colors automatically (Figure 1A). Becn1, as the biggest node and deepest color circle, was assigned the largest weight. Top 10 Hubba nodes were identified, which revealed the nearest interactions among these five genes, including Map1lc3b, Pik3r4, Atg14, Atg16l1, Atg7, Ambra1, Sqstm1, Atg9a, Atg5, and Atg13 (Figure 1B).

Figure 1 The construction of PPI network. According to the five coexpressed differential genes, other sixty interactive genes were predicted (A). The top 10 genes constructed a circle network using CytoHubba plug-in of Cytoscape software (B). The depth of the color and the size of the node correspond to the weighted score.

Visualization of Functional Enrichment and Interactive Networks

To elucidate the function of these 65 genes, GO annotations and KEGG pathway analysis were performed using the Metascape tool. GO functional enrichment results showed that DEGs were mainly participated in different types of autophagy. As shown in Figure 2, the enriched biological process (BP) terms were autophagy (GO:0006914), selective autophagy (GO:0061912), and lysosomal microautophagy (GO:0016237); The enriched cellular component (CC) terms were autophagosome (GO: 0005776), phosphatidylinositol 3-kinase complex, class III (GO:0035032), Bcl-2 family protein complex (GO:0097136), and Atg1/ULK1 kinase complex (GO:1990316); The enriched molecular function (MF) terms were ubiquitin protein ligase binding (GO:0031625), protein kinase binding (GO:0019901), BH3 domain binding (GO:0051434). KEGG pathway enrichment analysis revealed that these genes were involved in many signaling pathways, including autophagy-animal (ko04140), mitophagy-animal (ko04137), NOD-like receptor signaling pathway (ko04621), and longevity regulating pathway (mmu04211). The potential key GO and KEGG lists are referred to in Supplementary Table 2.

Figure 2 GO and KEGG enrichment analysis. A histogram of GO enrichment analysis of DEGs, included biological process (A), cellular component (B), and molecular function (C); the enrichment KEGG pathways (D).

To identify interactive hub genes, MCODE cluster analysis was used. Four key modules with different colors were enriched, such as MCODE1, MCODE2, MCODE3 and MCODE4 (Figure 3). MCODE1 mainly represented autophagy; MCODE2 included macroautophagy, and mitochondrion disassembly; MCODE3 included BH3-only proteins associate with and inactivate anti-apoptotic BCL-2 members and apoptosis; MCODE4 included positive regulation of cytokine production involved in immune response, and NF-kappa B signaling pathway.

Figure 3 Construction of interactive network. GO enrichment analysis was applied to each MCODE network. The same color nodes represent an interactive network and perform similar biological functions (A). Four MCODE components were constructed with the screened hub genes (B). Each interactive network with different colors own different score values.

Identification of LIRCPs and Prediction of Short Linear Sequence Motifs

Atg8-family proteins, as core autophagic proteins, are widely studied. In order to identify Atg8 interacting proteins and their binding linear sequences, the iLIR database was performed to identify novel putative LICRPs. We found these five DEGs could be as LIRCPs to interconnect with ATG8-family proteins. Each gene corresponded to validated or predicted short linear sequence motifs with different PSSM scores. The first two protein sequences of Becn1 had already been validated, including TSFKIL and NSFTLI sequences. The other short linear sequence motifs of DEGs were predicted and shown in Table 2. The potential overlapping sequences were summarized and marked in Supplementary Table 3. Overall, these five genes may regulate autophagy by interacting with Atg8-family proteins.

Table 2 Identification of LIR Motifs

Prediction of miRNA and Build miRNA-mRNA Interaction

To determine upstream miRNAs that could regulate mRNA, a miRNA-mRNA interaction network was constructed. In total, 1026 miRNA-mRNA pairs were identified and generated a complex interaction network with a filter value of 1 in the miRWalk database. Based on the analysis outcome, each miRNA could regulate the expression of multiple genes, such as mmu-miR-709 and mmu-miR-6391, which corresponded to target genes among Becn1, Atg3, Atg12 and Gabarapl2, except Pik3c3. In addition, each gene can be regulated by individual miRNA. The top 20 selected interactive miRNAs were shown in Figure 4. More detailed information about each miRNA is listed in Supplementary Table 4.

Figure 4 The prediction of miRNA-mRNA interactive networks. The miRWalk database was used to predict and constructed interactive networks (A). The top 20 interactive miRNAs were screened out and analyzed with the CytoHubba plug-in (B).

Key Genes Validation in vitro and in vivo

To induce the OA model in mice, DMM surgery was performed. All knee joint samples were collected at 4 weeks after surgery. We observed uneven cartilage surfaces and partial cartilage lesions by HE staining in the OA group (Figure 5A). Meanwhile, Safranin O-Fast Green staining revealed reduced proteoglycan and rough surfaces in articular cartilage (Figure 5B) with a higher OARSI score (Figure 5C), suggesting cartilage matrix degeneration in OA mice rather than the sham mice. In other words, the DMM surgery induced OA model could be applied at 4 weeks postoperation. To clarify the association between autophagy and activated osteoclasts in subchondral bone, TRAP staining were observed in two groups. We found that an increase of TRAP positive osteoclasts in the OA model compared to the sham group (Figure 5D).

Figure 5 Experiential verification in vitro and in vivo. HE staining (A), Safranin O-Fast Green staining (B) and OARSI scores (C) were shown in the sham group and the OA group. TRAP staining in tibial subchondral bone between two groups (D). Comparison of mRNA expression levels between the two groups (E). Immunohistochemical observation of these five DEGs between two groups (F). Immunofluorescence staining of Becn1 (G) and analysis of fluorescence intensity (H) in the OA group compared to the sham group. Asterisk (*) representing statistical difference is p<0.05.

To explore the role of autophagy in osteoclast formation between the sham and OA groups, BMMs were cultured and induced with M-CSF and RANKL in vitro. And the mRNA levels of these five key autophagy-related genes were examined using RT–PCR. As shown in Figure 5E, all five genes were upregulated.

To further validate these five differential genes expression in vivo, IHC staining and IF staining were used to evaluate. According to the IHC results, there shown the positive expression of Becn1, Atg3, Atg12, Pik3c3, and Gabarapl1 in the OA group (Figure 5F). Furthermore, IF staining revealed that Becn1 was significantly more highly expressed at site of the subchondral bone of the OA group than the sham group (Figure 5G and H). In summary, these findings demonstrate there is a high correlation between autophagy gene and osteoclast activation in a mouse OA model.

Discussion

Osteoarthritis (OA) is a degenerative whole joint disease that involves progressive loss of the articular cartilage and subchondral bone tissue remodeling, while the exact pathological mechanism is still unclear. Increasing evidence suggests that abnormal activation of osteoclasts in the subchondral bone mediates the progression of OA, and the inhibition of osteoclast activity could modify subchondral bone remodeling and attenuate articular cartilage degeneration.24 Previous studies pointed out that multiple autophagy-related genes and proteins had an impact on the differentiation and function of osteoclasts.25 And autophagy could increase the activity of osteoclasts and stimulate osteoclast-mediated bone resorption in vitro and in vivo.12 However, there are limited studies that have focused on the autophagy role of osteoclasts in the development of OA. Research on the autophagy regulatory network of osteoclasts would provide a new idea to elaborate the specified mechanism of OA. To mine and identify biomarkers and related pathways, the ways of bioinformatics analysis and experimental validation might open the door to study OA diseases.

Relying on the search of GEO databases, GSE41342 and GSE57468 were selected to analyze the underlying autophagy genes in osteoclasts by a bioinformatics approach. Five key genes were identified by GEO2R and ScanGEO analysis tools, such as Becn1, Atg3, Atg12, Pik3c3, and Gabarapl2, indicating that these 5 genes were as coexpressed autophagy-related genes of osteoclasts in OA disease. Among these genes, Becn1 is a key regulator of autophagy. Becn1, as a core component of the Pik3c3, induces the formation of autophagosomes in mammalian systems.26 And Becn1 has been shown to regulate osteoclast differentiation and then mediate bone homeostasis.12 Taken together, Becn1, as an important target, could be involved in every main step in the autophagic pathway; Atg3 is an ubiquitin-like–conjugating enzyme essential for autophagy. It plays a role in the regulation of autophagosome formation, nucleation and elongation of phagophores.27 However, it is still unknown how the Atg3 gene regulates osteoclasts; Atg12 is involved in the formation of autophagic vesicles, and it can be coupled to Atg5 through the ubiquitin-like coupling system to regulate the process of autophagy.28 A previous study suggested Atg12 and Atg5 promote osteoclast differentiation along with autophagic flux under hypoxia stimulation;29 Gabarapl2, a member of the Gabarap family, is involved in autophagosomal maturation and degradation.30 Especially in mitochondrial autophagy, Gabarapl2 could eliminate mitochondria to a basic level to meet cellular energy requirements and prevent excessive ROS production;31 Pik3c3, as a catalytic subunit of the Pik3 complex, are essential for autophagosome initiation and maturation.32 Pik3c3-c1 is a key regulator of autophagosome formation, and Pik3c3-c2 is required for maturation and sorting of autophagosomes and endosomes as they are trafficked to the lysosome. Although little information is available regarding Pik3c3 in osteoarthritis with autophagy pathways, it may be a novel therapeutic target.

These 5 core genes were imported and deduced to construct a PPI network. Sixty-five genes in total were acquired and the PPI suggested that the nearest signal transduction components were connected with the above 5 genes. We calculated that the largest weighted node was Becn1, and found Map1lc3b and Sqstm1 occupied relatively important positions in the top 10 genes. Some of these genes have been confirmed to be involved in the autophagy of osteoclasts and the pathophysiology of OA, such as Pi3k, Ulk1, and Bcl2l1.33 The other autophagy-related proteins predicted in PPI, such as Atg5, Atg7, and Atg4b, were essential for OC ruffled border generation, secretory activity and bone resorption.25 Atg5 is necessary for LC3-II to target the ruffled border and the presence of LC3-II in the membrane could promote fusion with secretory lysosomes which is required for bone resorption. Similar to Atg5, the lack of Atg7 suppressed the conversion of LC3I to LC3II and inhibited the absorption of osteoclasts. Knockdown of Atg7 inhibits the expression of TRAP and cathepsin K (CTSK), indicating that autophagy is involved in OC precursor differentiation.13 That is to say, the construction of a PPI network is feasible and meaningful, and these five genes may be targets for subsequent research.

According to the predicted PPI network of DEGs, the enrichment analysis results came from Metascape. Autophagy, selective autophagy, and lysosomal microautophagy were regarded as the most crucial biological processes. Autophagosome, Bcl-2 family protein complex, and Atg1/ULK1 kinase complex acted as cellular component. Regarding KEGG enrichment, autophagy-animal, mitophagy-animal, ferroptosis signaling pathway, NOD-like receptor signaling pathway, and NF-kappa B signaling pathway were deduced. Meanwhile, Metascape applies the mature recognition algorithm called MCODE to extract protein complexes and infer more biologically interpretable results. This analysis helps to identify more closely related genes in the network. The autophagy-other, macroautophagy and apoptosis pathways are closely related to different hub genes. Besides autophagy-related pathways, recent studies have shown that NOD-like receptor X1 (NLRX1) significantly impacts various diseases, including osteoarthritis.34 NOD-like receptors (NLRs) are a family of innate immune receptors that play key roles in host defense and inflammation. Stimulation of NOD2 leads to increased autophagosome formation by recruiting the autophagy protein ATG16L1 to the plasma membrane.35 How it recruits ATG16L1 and how it is transported to the plasma membrane are unknown. The Bcl-2 protein family, which contains at least one Bcl-2 homology (BH) region, plays a dual role in the regulation of autophagy.36 Under normal circumstances, Beclin1 binds to the antiapoptotic proteins Bcl-2 to inhibit autophagy, whereas under nutrient starvation, proapoptotic BH3-only proteins from the Bcl-2 family, such as Bad, can competitively disrupt the interaction between Beclin1 and Bcl-2/Bcl-XL to induce autophagy.37 Bcl-2, as a bridge between autophagy and apoptosis, has a large space to explore the activated osteoclast mechanism in OA. The relationship between ferroptosis and autophagy in abnormally activated osteoclasts is still unclear. Some hypotheses have suggested that the autophagic machinery can trigger ferroptosis by degradation ferritin38 or that ferroptosis-inducing conditions can stimulate autophagy.39 The crosstalk between ferroptosis and autophagy is worthy of further study, especially in OA disease.

Using the iLIR database, we identified the known protein sequences and predicted unknown LC3-interaction regions among these five proteins. These five proteins could regulate the process of autophagy associated with Atg8 family proteins. Atg8 proteins, as the prototype of the family, occupy a central position to regulate autophagic machinery.40 Atg8 proteins not only are essential for the elongation and closure of the phagophore, but also interact with selective autophagy-related receptor mediated autophagic degradation by targeting specific substrates.41 Another new mechanism suggested there exists an the interaction between receptors and Atg8-family proteins by an LC3-interacting region.42 These five genes were identified as LIR motif-containing proteins, and all the shortest sequences available in the database were mined. Potential binding sequences are disordered and influence the regulation of autophagy. The potential normal and disordered interaction protein sequences might provide precise targets to intervene.

miRNAs are short noncoding RNAs that regulate multiple cellular processes, including autophagy.43 To explore miRNA-mediated regulation of autophagy genes, miRNAs were predicted and a miRNA-mRNA interactive network was constructed. Each gene corresponded to individual miRNAs, and partial miRNAs were connected with different genes. Mmu-miR-709 and mmu-miR-6391 interacted with target genes among Becn1, Atg3, Atg12 and Gabarapl2. These results suggest that mmu-miR-709 and mmu-miR-6391 may play an important role in regulating autophagy in osteoclasts, although they have not been reported in the relevant literature. These significantly differentially expressed miRNAs had potential binding sites for autophagy-related genes.

To verify the expression of these autophagy-related genes in OA, the DMM surgery was performed to establish an OA model. Four weeks after surgery, osteoarthritis syndrome was observed in the OA group, and the reliability was verified by DMM surgery. As previously described, early changes in articular cartilage and subchondral bone of osteoarthritis occurred at 2 weeks postoperation, and the late stage osteoarthritic lesions at 6–10 weeks postsurgery.44–46 Here, we selected 4 weeks as the time point to observe that corresponded to the time of the GEO database, and the outcome suggested there that osteoarthritic lesions existed. Then, osteoclast precursor cells were induced to differentiate into osteoclasts in vitro, and the higher expression of these autophagy genes in the OA group was determined. These outcomes indicated autophagy in osteoclasts could be involved in the pathologies of OA at 4 weeks after surgery.

Many previous studies have focused on autophagy regulating chondrocyte function to ameliorate osteoarthritis,47 or autophagy regulates the differentiation of osteoclasts to reduce bone erosion with a rheumatoid arthritis model,48 few studies have shown the aberrant autophagy of osteoclasts has been linked to osteoarthritis. Research on the autophagy function of osteoclasts might provide an alternative idea to elaborate mechanics in OA.

Our study has certain limitations. To determine the molecular mechanisms of the pathogenesis of OA, an OA model from mice rather than clinical human tissues was selected for analysis. Since it is difficult to acquire OA samples and normal joint bone tissue during surgery at the early or middle stage, surgical OA models in mice that mimicked symptoms of posttraumatic OA in humans,49 were safe and effective, which induced rapid and severe joint degeneration after surgery.50 Furthermore, our study is also limited by observing autophagy of osteoclasts in subchondral bone at only one time point. More experimental validation with different time points after surgery should be performed to identify potential biomarkers and targets.

Conclusion

Bioinformatics analysis including identified key genes and pathways, then predicted iLIR motifs combined with mRNA-miRNA interactions, and experimental validation of differential genes in vitro and in vivo could be as approaches to explore the role of autophagy in osteoclasts. In summary, our research provided an important molecular mechanism from the perspective of autophagy in osteoclasts focused on OA disease, which might provide novel insight for OA research.

Abbreviations

OA, Osteoarthritis; GEO, Gene Expression Omnibus; DEGs, Differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, Biological process; CC, Cellular component; MF, Molecular function; PPI, Protein–protein interaction; LIR, LC3-interacting region; LIRCPs, LC3-interacting region containing proteins; PSSM, Position-specific scoring matrix; DMM, Destabilization of the medial meniscus; BMMs, Bone marrow-derived macrophage cells; HE, Hematoxylin and eosin; TRAP, Tartrate-resistant acid phosphatase; RT–PCR, Real-time polymerase chain reaction; IHC, Immunohistochemistry; IF, Immunofluorescence; OARSI, Osteoarthritis Research Society International; MCC, Maximal Clique Centrality; NLRX1, NOD like receptor X1; NLRs, NOD-like receptors; BH, Bcl-2 homology.

Ethics Approval and Consent to Participate

Ethical approval was obtained from the First Affiliated Hospital of Soochow University. All methods and experimental protocols were carried out in accordance with the guidelines and regulations of the Care and Use of Laboratory Animals of Soochow University.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: this work is supported by grants from the National Natural Science Foundation of China (Nos. 82072425, 82072498, 81873991, 81873990, 81672238 and 81472077), the Young Medical Talents of Jiangsu Province (No. QNRC2016751), the Natural Science Foundation of Jiangsu Province (Nos. BK20180001 and BE2021650), and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD) and Special Project of Diagnosis and Treatment Technology for Key Clinical Diseases in Suzhou (LCZX202003).

Disclosure

The authors declare that they have no conflicts of interest for this work.

References

1. Katz JN, Arant KR, Loeser RF. Diagnosis and treatment of hip and knee osteoarthritis: a review. JAMA. 2021;325(6):568–578. doi:10.1001/jama.2020.22171

2. Hunter DJ, Bierma-Zeinstra S. Osteoarthritis. Lancet. 2019;393(10182):1745–1759. doi:10.1016/S0140-6736(19)30417-9

3. Martel-Pelletier J, Barr AJ, Cicuttini FM, et al. Osteoarthritis. Nat Rev Dis Primers. 2016;2:16072. doi:10.1038/nrdp.2016.72

4. Teitelbaum SL. Bone resorption by osteoclasts. Science. 2000;289(5484):1504–1508. doi:10.1126/science.289.5484.1504

5. Schett G, Hayer S, Zwerina J, et al. Mechanisms of disease: the link between RANKL and arthritic bone disease. Nat Clin Pract Rheumatol. 2005;1(1):47–54. doi:10.1038/ncprheum0036

6. Jaiprakash A, Prasadam I, Feng JQ, et al. Phenotypic characterization of osteoarthritic osteocytes from the sclerotic zones: a possible pathological role in subchondral bone sclerosis. Int J Biol Sci. 2012;8(3):406–417. doi:10.7150/ijbs.4221

7. Zhen G, Cao X. Targeting TGFβ signaling in subchondral bone and articular cartilage homeostasis. Trends Pharmacol Sci. 2014;35(5):227–236. doi:10.1016/j.tips.2014.03.005

8. Zhu S, Zhu J, Zhen G, et al. Subchondral bone osteoclasts induce sensory innervation and osteoarthritis pain. J Clin Invest. 2019;129(3):1076–1093. doi:10.1172/JCI121561

9. Ahn SH, Chen Z, Lee J, et al. Inhibitory effects of 2N1HIA (2-(3-(2-Fluoro-4-Methoxyphenyl)-6-Oxo-1(6H)-Pyridazinyl)-N-1H-Indol-5-Ylacetamide) on osteoclast differentiation via suppressing cathepsin K expression. Molecules. 2018;23(12):3139. doi:10.3390/molecules23123139

10. Shapiro IM, Layfield R, Lotz M, et al. Boning up on autophagy: the role of autophagy in skeletal biology. Autophagy. 2014;10(1):7–19. doi:10.4161/auto.26679

11. Mizushima N, Levine B, Cuervo AM, et al. Autophagy fights disease through cellular self-digestion. Nature. 2008;451(7182):1069–1075. doi:10.1038/nature06639

12. Arai A, Kim S, Goldshteyn V, et al. Beclin1 modulates bone homeostasis by regulating osteoclast and chondrocyte differentiation. J Bone Miner Res. 2019;34(9):1753–1766. doi:10.1002/jbmr.3756

13. Pierrefite-Carle V, Santucci-Darmanin S, Breuil V, et al. Autophagy in bone: self-eating to stay in balance. Ageing Res Rev. 2015;24(Pt B):206–217. doi:10.1016/j.arr.2015.08.004

14. Zhang Y, Cui Y, Wang L, et al. Autophagy promotes osteoclast podosome disassembly and cell motility athrough the interaction of kindlin3 with LC3. Cell Signal. 2020;67:109505. doi:10.1016/j.cellsig.2019.109505

15. Lin NY, Beyer C, Giessl A, et al. Autophagy regulates TNFα-mediated joint destruction in experimental arthritis. Ann Rheum Dis. 2013;72(5):761–768. doi:10.1136/annrheumdis-2012-201671

16. Duan R, Xie H, Liu ZZ. The role of autophagy in osteoarthritis. Front Cell Dev Biol. 2020;8:608388. doi:10.3389/fcell.2020.608388

17. Koeppen K, Stanton BA, Hampton TH, Wren J. ScanGEO: parallel mining of high-throughput gene expression data. Bioinformatics. 2017;33(21):3500–3501. doi:10.1093/bioinformatics/btx452

18. Szklarczyk D, Morris JH, Cook H, et al. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 2017;45(D1):D362–d368. doi:10.1093/nar/gkw937

19. Saito R, Smoot ME, Ono K, et al. A travel guide to Cytoscape plugins. Nat Methods. 2012;9(11):1069–1076. doi:10.1038/nmeth.2212

20. Zhou Y, Zhou B, Pache L, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019;10(1):1523. doi:10.1038/s41467-019-09234-6

21. Jacomin AC, Samavedam S, Promponas V, et al. iLIR database: a web resource for LIR motif-containing proteins in eukaryotes. Autophagy. 2016;12(10):1945–1953. doi:10.1080/15548627.2016.1207016

22. Dweep H, Sticht C, Pandey P, et al. miRWalk–database: prediction of possible miRNA binding sites by “walking” the genes of three genomes. J Biomed Inform. 2011;44(5):839–847. doi:10.1016/j.jbi.2011.05.002

23. Glasson SS, Blanchet TJ, Morris EA. The surgical destabilization of the medial meniscus (DMM) model of osteoarthritis in the 129/SvEv mouse. Osteoarthr Cartil. 2007;15(9):1061–1069. doi:10.1016/j.joca.2007.03.006

24. Zhang RK, Li GW, Zeng C, et al. Mechanical stress contributes to osteoarthritis development through the activation of transforming growth factor beta 1 (TGF-β1). Bone Joint Res. 2018;7(11):587–594. doi:10.1302/2046-3758.711.BJR-2018-0057.R1

25. Yin X, Zhou C, Li J, et al. Autophagy in bone homeostasis and the onset of osteoporosis. Bone Res. 2019;7:28. doi:10.1038/s41413-019-0058-7

26. Mei Y, Glover K, Su M, et al. Conformational flexibility of BECN1: essential to its key role in autophagy and beyond. Protein Sci. 2016;25(10):1767–1785. doi:10.1002/pro.2984

27. Singh AK, Pandey RK, Shaha C, et al. MicroRNA expression profiling of Leishmania donovani -infected host cells uncovers the regulatory role of MIR30A-3p in host autophagy. Autophagy. 2016;12(10):1817–1831. doi:10.1080/15548627.2016.1203500

28. Hanada T, Noda NN, Satomi Y, et al. The Atg12-Atg5 conjugate has a novel E3-like activity for protein lipidation in autophagy. J Biol Chem. 2007;282(52):37298–37302. doi:10.1074/jbc.C700195200

29. Zhao Y, Chen G, Zhang W, et al. Autophagy regulates hypoxia-induced osteoclastogenesis through the HIF-1α/BNIP3 signaling pathway. J Cell Physiol. 2012;227(2):639–648. doi:10.1002/jcp.22768

30. Hervouet E, Claude-Taupin A, Gauthier T, et al. The autophagy GABARAPL1 gene is epigenetically regulated in breast cancer models. BMC Cancer. 2015;15:729. doi:10.1186/s12885-015-1761-4

31. Antón Z, Landajuela A, Hervás JH, et al. Human Atg8-cardiolipin interactions in mitophagy: specific properties of LC3B, GABARAPL2 and GABARAP. Autophagy. 2016;12(12):2386–2403. doi:10.1080/15548627.2016.1240856

32. Brier LW, Ge L, Stjepanovic G, Thelen AM, Hurley JH, Schekman R. Regulation of LC3 lipidation by the autophagy-specific class III phosphatidylinositol-3 kinase complex. Mol Biol Cell. 2019;30(9):1098–1107. doi:10.1091/mbc.E18-11-0743

33. Fu L, Wu W, Sun X, et al. Glucocorticoids enhanced osteoclast autophagy through the PI3K/Akt/mTOR signaling pathway. Calcif Tissue Int. 2020;107(1):60–71. doi:10.1007/s00223-020-00687-2

34. Nagai-Singer MA, Morrison HA, Allen IC, et al. NLRX1 is a multifaceted and enigmatic regulator of immune system function. Front Immunol. 2019;10:2419. doi:10.3389/fimmu.2019.02419

35. Travassos LH, Carneiro LA, Ramjeet M, et al. Nod1 and Nod2 direct autophagy by recruiting ATG16L1 to the plasma membrane at the site of bacterial entry. Nat Immunol. 2010;11(1):55–62. doi:10.1038/ni.1823

36. Yang Z, Klionsky DJ. Mammalian autophagy: core molecular machinery and signaling regulation. Curr Opin Cell Biol. 2010;22(2):124–131. doi:10.1016/j.ceb.2009.11.014

37. Lomonosova E, Chinnadurai G. BH3-only proteins in apoptosis and beyond: an overview. Oncogene. 2008;27 Suppl 1(Suppl1):S2–S19. doi:10.1038/onc.2009.39

38. Zhou Y, Shen Y, Chen C, et al. The crosstalk between autophagy and ferroptosis: what can we learn to target drug resistance in cancer? Cancer Biol Med. 2019;16(4):630–646. doi:10.20892/j.issn.2095-3941.2019.0158

39. Gao M, Monian P, Pan Q, et al. Ferroptosis is an autophagic cell death process. Cell Res. 2016;26(9):1021–1032. doi:10.1038/cr.2016.95

40. Kalvari I, Tsompanis S, Mulakkal NC, et al. iLIR: a web resource for prediction of Atg8-family interacting proteins. Autophagy. 2014;10(5):913–925. doi:10.4161/auto.28260

41. Shpilka T, Weidberg H, Pietrokovski S, et al. Atg8: an autophagy-related ubiquitin-like protein family. Genome Biol. 2011;12(7):226. doi:10.1186/gb-2011-12-7-226

42. Noda NN, Ohsumi Y, Inagaki F. Atg8-family interacting motif crucial for selective autophagy. FEBS Lett. 2010;584(7):1379–1385. doi:10.1016/j.febslet.2010.01.018

43. Wang S, Deng Z, Ma Y, et al. The role of autophagy and mitophagy in bone metabolic disorders. Int J Biol Sci. 2020;16(14):2675–2691. doi:10.7150/ijbs.46627

44. Kim S, Han S, Kim Y, et al. Tankyrase inhibition preserves osteoarthritic cartilage by coordinating cartilage matrix anabolism via effects on SOX9 PARylation. Nat Commun. 2019;10(1):4898. doi:10.1038/s41467-019-12910-2

45. Fang H, Huang L, Welch I, et al. Early changes of articular cartilage and subchondral bone in the DMM mouse model of osteoarthritis. Sci Rep. 2018;8(1):2855. doi:10.1038/s41598-018-21184-5

46. Che X, Chi L, Park CY, et al. A novel method to detect articular chondrocyte death during early stages of osteoarthritis using a non-invasive ApoPep-1 probe. Arthritis Res Ther. 2015;17:309. doi:10.1186/s13075-015-0832-x

47. Takahashi K, Nakamura H, Ozawa H, et al. Effectiveness of radiofrequency hyperthermia for treating cartilage in guinea pigs with primary osteoarthritis. Cartilage. 2018;9(1):71–79. doi:10.1177/1947603516678974

48. Wu DJ, Gu R, Sarin R, et al. Autophagy-linked FYVE containing protein WDFY3 interacts with TRAF6 and modulates RANKL-induced osteoclastogenesis. J Autoimmun. 2016;73:73–84. doi:10.1016/j.jaut.2016.06.004

49. Glasson SS. In vivo osteoarthritis target validation utilizing genetically-modified mice. Curr Drug Targets. 2007;8(2):367–376. doi:10.2174/138945007779940061

50. Bendele AM. Animal models of osteoarthritis. J Musculoskelet Neuronal Interact. 2001;1(4):363–376.

Creative Commons License © 2022 The Author(s). This work is published and licensed by Dove Medical Press Limited. The full terms of this license are available at https://www.dovepress.com/terms.php and incorporate the Creative Commons Attribution - Non Commercial (unported, v3.0) License. By accessing the work you hereby accept the Terms. Non-commercial uses of the work are permitted without any further permission from Dove Medical Press Limited, provided the work is properly attributed. For permission for commercial use of this work, please see paragraphs 4.2 and 5 of our Terms.