Two independent human knee articular cartilage tissue mRNA expression profile datasets were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE114007 RNA-seq dataset included 20 OA samples and 18 normal samples from two different platforms, GPL18573 and GPL11154 (18). GSE169077 microarray dataset (GPL96 platform) included six OA samples and five normal samples, which served as a validation dataset for verifying hub genes. Detailed information on the datasets is shown in Table S1 . In addition, 222 ARGs were obtained from the HADb database (http://www.autophagy.lu). The workflow of this research is shown in Figure 1 .

The GSE114007 gene expression of profile was quality assessed using the “factoextra” package (https://cloud.r-project.org/package=factoextra/) in R software before analyzing differentially expressed genes (DEGs) in OA, and sample clustering results suggested that there were large differences between two different platforms ( Supplementary Figure S1 ); therefore, only data from the GPL18573 platform, which included 10 OA samples and 10 normal samples, were taken for subsequent analysis. The counts’ data were standardized by the transcripts per kilobase per million mapped reads (TPM), and then the DEGs between OA and normal knee cartilage tissue were screened using the R package “DESeq2” with the threshold set as |log2FoldChange| > 1 and p-value <0.05 (19). Subsequently, differentially expressed ARGs (DEARGs) were obtained by intersecting genes between DEGs and ARGs. The Venn diagram drawn through the “VennDiagram” package displayed the number of DEARGs (https://CRAN.R-project.org/package=VennDiagram). Heatmap, volcano plot, and boxplot were drawn through the “pheatmap”, “ggplot2”, and “ggpubr” packages of R software, respectively.

The Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment of DEARGs were conducted in the R package “clusterProfiler” (20), and gene IDs were converted using Perl scripts. Molecular function (MF), biological process (BP), and cellular component (CC) constitute the GO annotation. Under the conditions of p < 0.05 and q < 0.05, GO terms and signaling pathways with significant differences were screened, and then the R software packages “enrichment plot”, “ggplot2”, and “GOplot” (21) were used to present the results, respectively.

The correlation between DEARG was investigated using Spearman’s correlation in the “Corrplot” package. The STRING online database (http://string-db.org/) was used to perform the protein–protein interaction (PPI) network analysis and investigate the underlying relationships between DEARGs with an interaction score >0.4, and then Cytoscape (version 3.8.1) was used to analyze and visualize the screened networks. MCODE and CytoHubba, two plug-ins of Cytoscape, use different algorithms to obtain autophagy-related hub genes. CytoHubba was used to calculate hub genes based on four different algorithms: maximal clique centrality (MCC), degree, closeness, and maximum neighborhood component (MNC). The top 10 genes were screened out of each algorithm as autophagy-related hub genes. Finally, the hub genes of DEARGs were defined as the overlapped genes obtained by five algorithms, and these genes were screened by drawing an UpSet diagram through the R package “UpSetR”.

Single-sample gene set enrichment analysis (ssGSEA) was used to analyze the abundance of immune cell infiltration in OA and normal cartilage tissues and identify 28 types of immune cell infiltration (22). Immune cell infiltration between OA and normal samples was compared using the Wilcoxon test, and the results were then visualized by violin diagrams drawn through the “ggplots” package. The correlation analysis of 28 kinds of infiltrating immune cells was performed, and then the results were visualized using the “ggcorrplot” package. Finally, Spearman’s correlation analysis was used to analyze the correlation between hub ARGs and the infiltration degree of different immune cells, and the results were visualized using the “ggcorrplot” and “ggstatsplot” packages. Furthermore, CIBERSORT algorithm was also adopted to investigate the immune cell infiltration.(As shown in Supplementary Figure S2 ).

To understand the expression of these hub ARGs in cartilage, validation was performed on the GSE169077 OA dataset using the “limma” package in R software. The heatmaps and boxplots were performed using the “pheatmap” and “ggplot2” packages in R language.

To confirm the hub ARGs, knee cartilage tissue samples from five cases of OA and five cases of traumatic amputation were collected after obtaining written informed consent from all patients. The study was approved by the Ethics Committee of the First Affiliated Hospital of Guangxi Medical University (Nanning, China) and complied with the tenets of the Declaration of Helsinki.

The fresh cartilage tissues were rinsed with sterile phosphate-buffered saline (PBS) three times and then cut into the size of approximately 0.3–0.5-mm³ pieces. The tissues were first treated with trypsin/ethylenediaminetetraacetic acid (EDTA) (Solarbio, Beijing, China) at a concentration of 0.25% for 30 min and then digested with collagenase II (Solarbio, Beijing, China) at a concentration of 2 mg/mL for 4 h at 37°C. Dulbecco’s Modified Eagle Medium (DMEM; Gibco, Shanghai, China) was used to culture chondrocytes after isolation, which consisted of 10% (v/v) fetal bovine serum (FBS; Tianhang, Zhejiang, China) and 1% (v/v) antibiotics (penicillin 10,000 U/ml and streptomycin 10,000 μg/ml, Solarbio, Beijing, China), cultured at an incubator with 5% CO₂ at 37°C. Passages 2–3 of chondrocytes were used for further experiments. Chondrocytes were divided into two groups: the control group (cultured with normal DMEM) and the treatment group (OA group). Chondrocytes treated with 10 ng/mL IL-1β (Solarbio, Beijing, China) for 24 hours were considered as the OA group. Then, total cellular RNA and protein were collected for further analysis.

According to the instructions of HiPure Total RNA Mini Kit (Magen, Guangzhou, China), total RNA was extracted from normal and OA group chondrocytes, and then reverse transcription with a Reverse Transcription kit (Takara, Dalian, China) was performed after checking the RNA concentration by a micro-spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The mRNA expression levels of the DNA damage-inducible transcript 3 (DDIT3), cyclin-dependent kinase inhibitor 1A (CDKN1A), fos proto-oncogene (FOS), vascular endothelial growth factor A (VEGFA), rela proto-oncogene (RELA), heat shock protein family A member 5 (HSPA5), microtubule-associated protein 1 light chain 3 beta (MAP1LC3B), MYC proto-oncogene (MYC), and heat shock protein family A member 8 (HSPA8) were analyzed by qRT-PCR. Table 1 shows the primer sequences for the main qRT-PCR used in this study. qRT-PCR was performed as previously reported (23).

According to the instructions of RIPA Lysis Buffer (Beyotime, Shanghai, China), total protein was extracted from normal and OA chondrocytes, and then the concentration of the extracted protein was determined by bicinchoninic acid (BCA) protein detection kit (Beyotime, Shanghai, China). The protein samples (40 μg/lane) were then separated by 10% sodium dodecyl sulfate (SDS)–polyacrylamide gels, and the blots were incubated overnight using antibodies against CDKN1A (1:1,000; Beyotime), DDIT3 (1:1,000; Beyotime), MAP1LC3B (1:1,000; Proteintech, Chicago, IL, USA), MYC (1:2,000; Proteintech), and GAPDH (1:5,000; Sangon Biotech, Shanghai, China) at 4°C. After washing with TBST, the membranes were incubated with goat anti-rabbit horseradish peroxidase (HRP)-conjugated secondary antibody (1:5,000; Sangon Biotech) at room temperature for 1 hour. The protein blots were detected with an enhanced chemiluminescence (ECL) system (Beyotime, Shanghai, China), and images were acquired with the Amersham Imager (Cytiva, Uppsala, Sweden). The intensity of these blots was quantified using ImageJ software (NIH, Bethesda, MD, USA).

Immunohistochemical (IHC) staining was conducted using a universal two-step detection kit (PV-9000; ZSGB-BIO, Beijing, China) following the manufacturer’s instructions. After being dewaxed, the knee articular cartilage tissue slices were placed in EDTA antigen repair solution at 95°C for 15 min for antigen repair. Endogenous peroxidase’s activity was blocked with peroxidase blocking reagent for 10 min at room temperature (RT). An appropriate amount of rabbit anti-CDKN1A (1:200; Beyotime), rabbit anti-DDIT3 (1:200; Beyotime), rabbit anti-MAP1LC3B (1:200, Proteintech), and rabbit anti-MYC (1:200, Proteintech) primary antibody working solution was added dropwise on the tissue slices and then incubated in a wet box overnight at 4°C. After rinsing with PBS, the appropriate amount of goat anti-mouse/rabbit IgG polymer labeled with enhanced enzyme was added dropwise on tissue slices and incubated at 37°C for 20 min. Subsequently, the tissue slices were visualized with DAB staining, and then counterstaining was achieved with hematoxylin.

Data are expressed as mean ± standard deviation (SD). Gene expression or protein levels in the two groups of chondrocytes were compared by unpaired Student’s t-test or Wilcoxon test. R software (version 4.1.3) and GraphPad Prism 8.0 (San Diego, CA, USA) were used for statistical analysis and plots. All correlation analyses were performed using Spearman’s method. Statistical significance required p < 0.05 (two-sided).

Through differential expression analysis of GPL18573 platform-derived data in the GSE114007 dataset, 2,567 genes were significantly expressed in OA compared with control cartilage tissue based on |log2Fold Change| > 1 with p < 0.05. Twenty-nine DEARGs were obtained from 222 ARGs, which intersected with 2,567 DEGs, and five upregulated genes and 24 downregulated genes were screened ( Figure 2A , Table 2 ). Volcano plot and heatmap were used to visualize the expressions of the 29 DEARGs ( Figures 2B, C ). Furthermore, the boxplot showed 29 DEARG expression patterns in OA compared to control cartilage tissues ( Figure 3 ).

GO and KEGG pathway enrichment analyses were performed on these DEARGs using R software to explore their potential biological functions and pathways. The results of GO functional analysis revealed that the most significant items of GO enrichment included response to nutrient levels, cellular response to external stimulus, cellular response to extracellular stimulus, and regulation of autophagy (biological process); autophagosome, mitochondrial outer membrane, organelle outer membrane, and outer membrane (cellular component); chaperone binding, ubiquitin-protein ligase binding, ubiquitin-like protein ligase binding, and heat shock protein binding (molecular function) ( Figures 4A, B ). The results of KEGG pathway enrichment analysis showed that DEARGs were mainly enriched in autophagy, mitophagy, and inflammation-related pathways ( Figures 4C, D ).

Spearman’s correlation analysis was used to explore the correlation of these DEARGs in OA. The results of the study indicated an interaction between 29 DEARGs ( Figure 5A ). A PPI network analysis was performed to determine the interactions between these DEARGs. PPI networks showed that DEARGs were involved in 26 nodes and 77 edges ( Figure 5B ) and the top 20 interactive genes ( Figure 5C ). According to the MCODE plug-ins of Cytoscape software, two clusters can be obtained ( Figure 5D ). The first cluster consists of 11 genes (FOS, CDKN1A, VEGFA, RELA, DDIT3, PPP1R15A, ERO1L, CDKN1B, DNAJB9, HSPA8, and CCL2), whereas the second cluster involves six genes (BNIP3, MYC, HSPA5, ULK1, BNIP3L, and MAP1LC3B).

Subsequently, the CytoHubba plug-in was used to calculate the PPI network constructed by the 26 DEARGs using four algorithms (MCC, MNC, Degree, and Closeness) and selected the top 10 genes as hub genes ( Figures 6A–D ). Finally, hub genes obtained by the CytoHubba algorithm were intersected with those calculated using the MCODE algorithm, and nine common hub genes were obtained, which included CDKN1A, DDIT3, FOS, VEGFA, RELA, MAP1LC3B, MYC, HSPA5, and HSPA8 ( Figure 6E ).

The expression pattern of hub genes was validated based on the GSE169077 dataset. The heatmap showed the expression levels of nine hub genes ( Figure 7A ). Additionally, as displayed in Figures 7B–J , the expression levels of CDKN1A, DDIT3, VEGFA, RELA, MAP1LC3B, MYC, and HSPA5 in OA cartilage were significantly lower than those in control samples, whereas HSPA8 was significantly upregulated (p < 0.05), which was in agreement with bioinformatics analysis of GSE114007 dataset. However, no significant difference was found in FOS expression.

ssGSEA was applied to explore immune cell infiltration in the OA cartilage. The results are shown in Figure 8A ; the proportions of regulatory T cells, activated dendritic cells, central memory CD4 T cells, T follicular helper cell, type 2 helper cell, macrophage, central memory CD8 T cells, gamma delta T cells, and immature dendritic cell were higher in the OA group and lower in the control group. In contrast, the proportions of activated CD56 bright nature killer cells, mast cells, type 17 T helper cells, B cells, eosinophil, and effector memory CD8 T cells were relatively low in the OA group compared to the control group. Subsequently, the correlation of 28 subpopulations of immune cell infiltrate was further analyzed. As shown in Figure 8B , except for immature B cells, a significant correlation was found among infiltrating immune cells. Finally, to explore the relationship between hub DEARGs and infiltrating immune cells, Spearman’s correlation analysis was performed ( Figure 8C ). The results suggest that the most positively correlated autophagy–immunocyte gene pair is FOS–Eosinophil, while the most negatively correlated autophagy–immunocyte gene pair is DDIT3–Immature dendritic cell ( Figures 8D, E ).

To verify the expression of hub DEARGs in OA chondrocytes, we further performed qRT-PCR and Western blotting. qRT-PCR results showed that the mRNA expressions of CDKN1A, DDIT3, FOS, MAP1LC3B, MYC, and HSPA8 were significantly low in the OA group; HSPA5 and VEGFA were significantly highly expressed (p < 0.05); and RELA gene expression was not significantly changed compared to normal controls ( Figure 9A ). Moreover, Western blotting results showed that the expression levels of CDKN1A, DDIT3, MAP1LC3B, and MYC were significantly decreased in the OA group ( Figure 9B ). Further, the IHC staining showed that compared with the normal group, expressions of CDKN1A, DDIT3, MAP1LC3B, and MYC were low in the OA articular cartilage tissue ( Figure 9C ).