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

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.

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 (), 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 (). 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 () 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 () 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 () 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 () 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 () 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 (), 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

Becn1	Atg3	Atg12	Pik3c3	Gabarapl2
Ulk3	Atg10	Atg5	Gabarapl1	Ifna11
Atg7	Atg16l1	Ulk1	Atg4d	Prkaa2
Prkaa1	Atg4c	Ifng	Atg16l2	Ifna9
Atg14	Atg4a	Ifnab	Ifna1	Ins1
Atg4b	Ins2	Ulk2	Ifna4	Pik3r4
Gabarap	Atg13	Ifna2	Ifna5

Note: A total of 34 DEGs were filtered out and the five genes in the first row were as coexpressed differential genes.

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 (). 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 .

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 . Overall, these five genes may regulate autophagy by interacting with Atg8-family proteins.

Table 2

Identification of LIR Motifs

Gene Name	Motif	Start	End	Pattern	PSSM Score
Becn1	WxxL	28	33	TSFKIL	8
	WxxL	93	98	NSFTLI	9
	WxxL	252	257	VRYAQI	2
	WxxL	278	283	GQFGTI	5
	WxxL	296	301	VEWNEI	14
	WxxL	324	329	QRYRLV	3
Atg3	WxxL	231	236	HMYEDI	10
	WxxL	287	292	HMYLLI	9
Atg12	WxxL	90	95	KKFLKL	8
	WxxL	100	105	QLFIYV	3
Pik3c3	xLIR	214	219	FMYLMV	13
	WxxL	5	10	EKFHYI	5
	WxxL	33	38	KSYKAV	8
	WxxL	44	49	LKFSGL	3
	WxxL	88	93	NEWLKL	17
	WxxL	95	100	VKYPDL	3
	WxxL	144	149	KVWPNV	6
	WxxL	196	201	LTFREI	7
	WxxL	212	217	SNFMYL	7
	WxxL	219	224	VEFRCV	3
	WxxL	229	234	KEYGIV	8
	WxxL	248	253	TSFELV	9
	WxxL	406	411	ENFDDI	11
	WxxL	492	497	YLYWYV	−1
	WxxL	515	520	EMYLNV	9
	WxxL	629	634	GKYPVI	6
	WxxL	668	673	TPYKVL	6
	WxxL	730	735	AGYCVI	8
	WxxL	760	765	IDFGYI	7
	WxxL	804	809	TAFLHL	6
	WxxL	810	815	RRYSNL	5
	WxxL	869	874	ALFAAV	2
Gabarapl2	WxxL	58	63	AQFMWI	3
	WxxL	75	80	AIFLFV	6
	WxxL	102	107	DGFLYV	11

Notes: These five genes predicted LC3-interacting regions with different positions and PSSM scores. Dark gray have already been validated and light grey represented not been confirmed.

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 .

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.