Identification of autophagy-related genes in osteoarthritis articular cartilage and their roles in immune infiltration

Background: Autophagy plays a critical role in the progression of osteoarthritis (OA), mainly by regulating inflammatory and immune responses. However, the underlying mechanisms remain unclear. This study aimed to investigate the potential relevance of autophagy-related genes (ARGs) associated with infiltrating immune cells in OA.

Methods: GSE114007, GSE169077, and ARGs were obtained from the Gene Expression Omnibus (GEO) database and the Human Autophagy database. R software was used to identify the differentially expressed autophagy-related genes (DEARGs) in OA. Functional enrichment and protein–protein interaction (PPI) analyses were performed to explore the role of DEARGs in OA cartilage, and then Cytoscape was utilized to screen hub ARGs. Single-sample gene set enrichment analysis (ssGSEA) was used to conduct immune infiltration analysis and evaluate the potential correlation of key ARGs and immune cell infiltration. Then, the expression levels of hub ARGs in OA were further verified by the GSE169077 and qRT-PCR. Finally, Western blotting and immunohistochemistry were used to validate the final hub ARGs.

Results: A total of 24 downregulated genes and five upregulated genes were identified, and these genes were enriched in autophagy, mitophagy, and inflammation-related pathways. The intersection results identified nine hub genes, namely, CDKN1A, DDIT3, FOS, VEGFA, RELA, MAP1LC3B, MYC, HSPA5, and HSPA8. GSE169077 and qRT-PCR validation results showed that only four genes, CDKN1A, DDT3, MAP1LC3B, and MYC, were consistent with the bioinformatics analysis results. Western blotting and immunohistochemical (IHC) showed that the expression of these four genes was significantly downregulated in the OA group, which is consistent with the qPCR results. Immune infiltration correlation analysis indicated that DDIT3 was negatively correlated with immature dendritic cells in OA, and FOS was positively correlated with eosinophils.

Conclusion: CDKN1A, DDIT3, MAP1LC3B, and MYC were identified as ARGs that were closely associated with immune infiltration in OA cartilage. Among them, DDIT3 showed a strong negative correlation with immature dendritic cells. This study found that the interaction between ARGs and immune cell infiltration may play a crucial role in the pathogenesis of OA; however, the specific interaction mechanism needs further research to be clarified. This study provides new insights to further understand the molecular mechanisms of immunity involved in the process of OA by autophagy.

1 Introduction

Osteoarthritis (OA) is a common chronic inflammatory disease in clinical practice, and its development is related to many factors including age, joint damage, and obesity (1). The development is characterized by inflammation, cartilage degeneration, narrowing of the joint space, formation of osteophytes, and sclerosis of the subchondral bone (2). Currently, there is no effective treatment to alleviate the disease. While non-steroidal anti-inflammatory drug (NSAID) interventions mainly address pain and inflammation, OA cannot be prevented. Joint replacement surgery is often a solution for end-stage patients but imposes a significant financial burden (3, 4). Therefore, understanding the molecular mechanisms of OA pathogenesis will provide new ideas for the treatment of OA. The process of OA has been reported to be associated with aging, inflammation, apoptosis, and autophagy (1, 5, 6). Previous studies have shown that the dysfunction of autophagy is a major factor in OA (7). Therefore, identifying novel molecular biological targets is crucial for the in-depth study of the underlying molecular mechanisms by which autophagy regulates the pathogenesis of OA.

Autophagy is an intracellular biological degradation system with highly evolutionary conservative features that slow down the OA process by restoring cellular dysfunction in cartilage (8). Autophagy has great clinical prospects as a potential therapeutic target for OA; however, the specific pathogenesis of OA exacerbated by dysregulated autophagy has not been elucidated. Previous research has indicated that mTOR, LC3-II, Beclin-1, and p62 were associated with OA progression and provide promising therapeutic targets for OA (9, 10). As an important key regulator, mTOR affects autophagic activity and is a key target for the autophagic pathway (8). Microtubule-associated proteins 1A/1B light chain 3B (LC3) are involved in the biosynthesis of autophagosomes. In OA, increased chondrocyte apoptosis was associated with the decrease of LC3 (11). The levels of autophagy-related markers, such as ATG5, LC3-II, and Beclin-1 gene expression, were significantly increased in chondrocytes in the early stage of OA. With the aggravation of the disease, oxidative stress-induced damage increased gradually while decreasing autophagy, resulting in chondrocyte hypertrophy (12). These research data demonstrate that autophagy plays a key role in OA. However, autophagy-related genes (ARGs) in OA development are still largely unknown. Therefore, in-depth exploration of autophagy-related markers in OA through bioinformatics could help identify new potential biomarkers for the treatment of osteoarthritis. In addition, studies have indicated that immune regulation may also play a key role in the pathogenesis of OA (13, 14). Autophagy is an important process in regulating immune responses, which not only eliminates infectious agents and modulates inflammatory responses but also selects antigen presentation and regulates T-cell homeostasis and activation (15, 16). Huang et al. showed that autophagy and immune regulation played different roles in OA and RA, but both are significantly correlated (17). However, few reports have comprehensively investigated the relationship between ARG expression and cartilage immune infiltration in OA and the role of ARGs’ biological functions in cartilage immune infiltration. Therefore, we hypothesized that ARGs were differentially expressed in OA and normal cartilage tissues, and the change of autophagy level may related to immune cell infiltration.

In this study, we identified ARGs associated with OA based on OA-related data in the Gene Expression Omnibus (GEO) database. We found the association of the hub genes including CDKN1A, DDIT3, MAP1LC3B, and MYC with immune cell infiltration in OA, which was further verified using qRT-PCR, Western blotting, and immunohistochemical staining. Our study hopes to identify ARGs in OA cartilage tissue and then determine whether they regulate OA by modulating immune cell infiltration, which will provide a novel perspective for OA pathogenesis.

2 Materials and methods

2.1 Acquisition of raw data and pre-processing

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.

FIGURE 1

Figure 1 The workflow chart of this study.

2.2 Identification of differentially expressed autophagy-related genes

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.

2.3 Functional enrichment analysis of DEARGs

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.

2.4 Correlation analysis and PPI network analysis of DEARGs

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”.

2.5 Immune cell infiltration and correlation analysis

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

2.6 Validation of hub gene expression with other osteoarthritis datasets

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.

2.7 Clinical specimen collection and chondrocyte isolation and culture

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.

2.8 RNA extraction and qRT-PCR

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

TABLE 1

Table 1 Primer sequences used in the qRT-PCR experiments.

2.9 Protein extraction and Western blotting

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

2.10 Immunohistochemical analysis

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.

2.11 Statistical analysis

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

3 Results

3.1 Differential expression of ARGs in osteoarthritis cartilage

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

FIGURE 2

Figure 2 Identification of DEARGs. (A) Veen diagram showing the DEARGs. (B) Volcano plot displaying significantly differentially expressed ARGs. Red dots represent the upregulated genes, and dark purple dots denote the downregulated genes, with thresholds of |log2FoldChange| ≥ 1 and adjusted p-value <0.05. (C) The expressions of 29 DEARGs in OA are displayed in the heatmap. Red signifies significantly upregulated ARGs, and purple indicates significantly downregulated ARGs in the samples. DEARGs, differentially expressed autophagy-related genes; ARGs, autophagy-related genes; OA, osteoarthritis.

TABLE 2

Table 2 The 29 DEARGs in OA cartilage tissues compared to healthy control.

FIGURE 3

Figure 3 Boxplot displaying differential expressed of 29 autophagy-related genes in OA and normal cartilage samples. p-Values were calculated using Wilcoxon test. *p < 0.05; **p < 0.01; ***p < 0.01. ns, not significant; OA, osteoarthritis.

3.2 GO and KEGG enrichment analyses of DEARGs

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

FIGURE 4

Figure 4 Functional enrichment analysis of differentially expressed autophagy-related genes. Bubble plot (A) and Circor chart (B) of Gene Ontology (GO) enrichment analysis results of 29 DEARGs in biological process (BP), cellular component (CC), and molecular function (MF). Bubble plot (C) and Circor chart (D) of KEGG pathway enrichment analysis by DEARGs. DEARGs, differentially expressed autophagy-related genes; KEGG, Kyoto Encyclopedia of Genes and Genomes.

3.3 Correlation analysis and PPI network analysis of DEARGs

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

FIGURE 5

Figure 5 Correction analysis and protein–protein interaction (PPI) analysis of the 29 differentially expressed autophagy-related genes. (A) Spearman’s correction analysis of the 29 differentially expressed autophagy-related genes. (B) PPI network of DEARGs; three disconnected nodes were removed. Circle size represents the node degree. (C) The interaction number of the top 20 DEARGs. (D) Two DEARG clusters were obtained using the MCODE plug-ins. DEARGs, differentially expressed autophagy-related genes.

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

FIGURE 6

Figure 6 Identification of hub genes. (A–D) The PPI network constructed by 26 DEARGs was calculated using four algorithms, namely, MCC, MNC, Degree, and Closeness, and the top 10 genes were selected as hub genes. (E) UpSet diagram displays the intersection of five different algorithms, resulting in nine shared hub genes. PPI, protein–protein interaction; DEARGs, differentially expressed autophagy-related genes.

3.4 Validation of the hub DEARGs by microarray dataset

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.

FIGURE 7

Figure 7 Cross-validation of hub genes. (A) The expression levels of CDKN1A, DDIT3, FOS, VEGFA, RELA, MAP1LC3B, MYC, HSPA5, and HSPA8 were verified by GSE169077 dataset, the results of which are presented as heatmap. (B–J) Detailed expression of nine hub genes. **p < 0.001. ns, not significant.

3.5 Assessment of immune cell infiltration and correlation analysis of hub DEARGs and infiltrating immune cells

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

FIGURE 8

Figure 8 Evaluation of immune cell infiltration and correlation analysis between hub genes and immune cell infiltration. (A) Violin plot of differential infiltrating fractions of all 28 immune cells between OA and normal samples. (B) Spearman’s correlation heatmap demonstrates the correlation of all 28 immune cells. The blank squares in the upper left corner represent p-values >0.05, red presents a positive correlation, and blue denotes a negative correlation. (C) Correlation heatmap displaying the correlations between nine hub genes and 28 infiltrating immune cells. The red and purple squares indicate that the hub genes have a significant correlation with infiltrating immune cells, and the blank squares represent p-values >0.05 for the association of hub genes with infiltrating immune cells. (D) The most positively correlated autophagy–immunocyte gene pair is FOS–Eosinophil. (E) The most negatively correlated autophagy–immunocyte gene pair is DDIT3–Immature dendritic cell. p-Values were calculated using the Wilcoxon test. ***p < 0.001; **p < 0.01; *p < 0.05. ns, not significant.

3.6 Validation of hub DEARGs using an in vitro OA cell model

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

FIGURE 9

Figure 9 Validation of hub DEARGs using an in vitro OA cell model. (A) Boxplots show mRNA expression levels of nine hub genes that were measured by qRT-PCR in OA and normal chondrocytes. (B) The expression levels of CDKN1A, DDIT3, MAPLC3B, and MYC were detected by Western blotting. (C) The expression levels of CDKN1A, DDIT3, MAPLC3B, and MYC were detected by IHC. *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001. ns, not significant; DEARGs, differentially expressed autophagy-related genes; OA, osteoarthritis; IHC, immunohistochemistry. Scale bar, 50 μm.

4 Discussion

OA is the most common chronic degenerative joint disease with limited treatment options, as its etiology has not been fully elucidated (4). Many studies in recent years have demonstrated that autophagy defects contribute to OA and age-related diseases; in contrast, autophagy activation could promote cell survival and reduce the severity of experimental OA (24, 25). However, the relationship between autophagy and OA is currently not fully understood. In addition, ARGs in OA cartilage have not been well explored using bioinformatics analysis. Therefore, we conducted this study to explore whether the expression of ARGs in OA cartilage tissue was significantly different from that in healthy controls. Additionally, we investigated the correlation between ARGs and immune infiltration, which can help understand the potential immune mechanisms of ARGs in cartilage and identify potential molecular targets in OA therapy.

In this study, bioinformatics analysis tools were used to identify ARGs and 29 potential DEARGs associated with OA. Among them, five genes were found to be upregulated in expression and 24 genes were downregulated in expression, as shown in Table 2. After construction using the PPI network, nine key ARGs (CDKN1A, DDIT3, FOS, VEGFA, RELA, MAP1LC3B, MYC, HSPA5, and HSPA8) were differentially expressed in OA cartilage. To understand the potential molecular biological function of DEARGs between OA and normal cartilage, we performed enrichment analysis using GO and KEGG. The enrichment results suggested that these DEARGs play important roles in various biological functions, including autophagy, mitochondrial autophagy, and inflammation-related pathways. In this case, we collected clinical samples to validate the expression of these nine genes. In vitro, cellular experiments showed that the expression of four genes, CDKN1A, DDIT3, MAP1LC3B, and MYC, were significantly consistent (p < 0.05) in bioinformatics analysis. Some ARGs have been previously studied in OA. CDKN1A is also known as protein p21. It was shown that CDKN1A expression was decreased in OA by inhibiting chondrocyte proliferation (26). Moreover, Huang et al. revealed that in cardiomyocytes, autophagy was inhibited when the expression level of p21 was decreased and that p21 regulates autophagy by interacting with MAP1LC3B(LC3B) (27). Yang et al. showed that DDIT3 knockdown suppressed the autophagy of chondrocytes in vitro and in vivo, causing a decrease in LC3B and Beclin1 gene expression; however, overexpression of DDIT3 significantly promoted autophagy in vitro (28). In microglia, DEX promotes NLRP3 inflammasome degradation via the autophagy-ubiquitin pathway and reduces MAP1LC3B expression, thereby reducing hippocampal inflammation (29). Anti-Dlx5 slows the progression of OA by downregulating the chondrocyte apoptosis-related gene MYC (30). These results are consistent with our bioinformatics analysis and experimental results. Moreover, four genes, namely, CDKN1A, DDIT3, MAP1LC3B, and MYC, may determine the progression of OA. The specific mechanisms need to be further investigated.

Many publications have reported that the infiltration of immune cells is a key factor in promoting the development of OA (31). For example, immune cells such as neutrophils, M1 macrophages, CD4+ T cells, and mast cells have a significant infiltration in OA synovium, suggesting that immune infiltration is a key target for the treatment of osteoarthritis (31–33). Many immune infiltration studies have been performed in OA synovial tissue, and relatively few immune infiltration studies have been performed in cartilage tissue. Therefore, we comprehensively evaluated the types of immune infiltration cells in OA cartilage tissues using ssGSEA, hoping to explore in-depth the potential mechanisms of immune infiltration in OA cartilage. We performed a correlation analysis of immune cell subsets in normal cartilage tissue and OA samples. It was found that the proportions of the 28 immune cell subpopulations differed significantly between healthy cartilage and OA cartilage samples. The proportion of regulatory T cells, activated dendritic cells, central memory CD4 T cells, T follicular helper cells, type 2 helper cells, macrophage, central memory CD8 T cells, gamma delta T cells, and immature dendritic cells were present in a higher proportion of OA samples compared with healthy controls. In contrast, the fraction of activated CD56 bright nature killer cells, mast cells, type 17 T helper cells, B cells, eosinophils, and effector memory CD8 T cells was lower in OA. Han et al. performed analysis using ssGSEA, MCPcounter, and ESTIMATE software and found that immune scores were significantly lower in normal control samples than in OA samples (34). The discrepancy could be caused by the difference in the methods of analysis and sample numbers. Previous studies have illustrated that the immune system is associated with autophagic pathways, including innate and adaptive immunity, and when autophagy is disturbed, it contributes to the implications of inflammatory diseases (35–37). However, few studies have shown the relationships between ARGs and immune cell infiltration in OA cartilage.

This study also analyzed the potential link between nine key ARGs and infiltrating immune cell subpopulations, contributing to an in-depth understanding of the regulatory role of the ARGs in OA cartilage and also providing new insights into the pathogenesis of OA. Our study showed that DDIT3 and FOS genes were significantly correlated with immune cells, with FOS expression being positively associated with eosinophils and DDIT3 expression being positively associated with immature dendritic cells. It has been reported that FOS may treat dermatomyositis by regulating the infiltration of immune cells (38). Deng et al. found that infiltration of eosinophils may be associated with OA progression (32). Dendritic cells (DCs) present a prominent immunomodulatory capacity in innate and adaptive immune responses and are specialized antigen-presenting cells. DCs in OA patients have been reported to secrete large amounts of inflammatory cytokines that exacerbate the inflammatory response (39). However, no report proved the role of eosinophils in OA; in addition, the pairs FOS–Eosinophil and DDIT3–Immature dendritic cell have also not been investigated in OA. These findings suggest that there may be an interaction between autophagy and the immune response in OA. All of these studies suggest that autophagy and immune infiltration may have a reciprocal regulatory role in OA cartilage. However, more experimental studies may be needed to elucidate in-depth the mechanisms of these ARG–immune cell interactions for validation.

However, this study still has some inevitable limitations. First, we re-mined and analyzed previously published datasets. Second, not multiple datasets but only one dataset was used for cross-validation, and datasets and clinical samples were insufficient in our study, which may lead to deviations in the results. Thus, larger sample sizes are required to confirm our findings. Finally, we simply verified the expression levels of hub DEARGs by qRT-PCR, Western blotting, and IHC staining in vitro by collecting clinical samples without in-depth exploration of the underlying mechanism. Therefore, further studies need to be conducted to reveal the underlying molecular mechanisms.

5 Conclusions

In conclusion, we identified nine key ARGs, namely, CDKN1A, DDIT3, FOS, VEGFA, RELA, MAP1LC3B, MYC, HSPA5, and HSPA8, in OA cartilage tissue by informatics analysis. Verified results confirm that CDKN1A, DDIT3, MAP1LC3B, and MYC could serve as potential biomarkers for OA cartilage tissue. Additionally, the correlation analysis between hub ARGs and immune infiltration in OA cartilage suggested that the interaction between ARGs and immune cell infiltration may be involved in regulating the pathogenesis of OA. Among them, the expression level of DDIT3 showed a strong negative correlation with immature dendritic cells. The present study revealed that autophagy may regulate the progression of OA by regulating immune infiltration, providing new insights into the molecular immune mechanisms of treating OA.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Ethics statement

The studies involving humans were approved by The Ethics Committee of Guangxi Medical University (NO. 2019-SB-058, 7 Mar. 2019). The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation in this study was provided by the participants’ legal guardians/next of kin.

Author contributions

JQ: Conceptualization, Data curation, Formal Analysis, Methodology, Resources, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. JZ: Conceptualization, Data curation, Methodology, Resources, Validation, Writing – original draft. J-JW: Conceptualization, Data curation, Software, Validation, Visualization, Writing – review & editing. XR: Data curation, Investigation, Writing – original draft. Q-LZ: Investigation, Writing – original draft. J-MZ: Funding acquisition, Project administration, Supervision, Writing – review & editing. N-HL: Funding acquisition, Project administration, Supervision, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was funded by the Guangxi Science and Technology Base and Talent Special Project, grant number GuikeAD19254003; the Guangxi Natural Science Foundation, grant number 2023GXNSFBA026032; the National Natural Science Foundation of China, grant number 81972120; and the Youth Science Foundation of Guangxi Medical University, grant number GXMUYSF202111.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1263988/full#supplementary-material

References

1. Loeser RF, Collins JA, Diekman BO. Ageing and the pathogenesis of osteoarthritis. Nat Rev Rheumatol (2016) 12(7):412–20. doi: 10.1038/nrrheum.2016.65

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Choi WS, Lee G, Song WH, Koh JT, Yang J, Kwak JS, et al. The ch25h-cyp7b1-roralpha axis of cholesterol metabolism regulates osteoarthritis. Nature (2019) 566(7743):254–8. doi: 10.1038/s41586-019-0920-1

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Gregori D, Giacovelli G, Minto C, Barbetta B, Gualtieri F, Azzolina D, et al. Association of pharmacological treatments with long-term pain control in patients with knee osteoarthritis: A systematic review and meta-analysis. JAMA (2018) 320(24):2564–79. doi: 10.1001/jama.2018.19319

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Hunter DJ, Bierma-Zeinstra S. Osteoarthritis. Lancet (2019) 393(10182):1745–59. doi: 10.1016/s0140-6736(19)30417-9

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Robinson WH, Lepus CM, Wang Q, Raghu H, Mao R, Lindstrom TM, et al. Low-grade inflammation as a key mediator of the pathogenesis of osteoarthritis. Nat Rev Rheumatol (2016) 12(10):580–92. doi: 10.1038/nrrheum.2016.136

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Castrogiovanni P, Ravalli S, Musumeci G. Apoptosis and autophagy in the pathogenesis of osteoarthritis. J Invest Surg (2020) 33(9):874–5. doi: 10.1080/08941939.2019.1576811

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Giuliani C. The flavonoid quercetin induces ap-1 activation in frtl-5 thyroid cells. Antioxidants (Basel) (2019) 8(5):112. doi: 10.3390/antiox8050112

PubMed Abstract | CrossRef Full Text | Google Scholar

8. 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

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Liu J, Fu Q, Liu S. Transcriptional regulation based on network of autophagy identifies key genes and potential mechanisms in human osteoarthritis. Cartilage (2021) 13(2_suppl):1431S–41S. doi: 10.1177/1947603520951632

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Meng Z, Shen B, Gu Y, Wu Z, Yao J, Bian Y, et al. Diazoxide ameliorates severity of experimental osteoarthritis by activating autophagy via modulation of the osteoarthritis-related biomarkers. J Cell Biochem (2018) 119(11):8922–36. doi: 10.1002/jcb.27145

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Carames B, Taniguchi N, Otsuki S, Blanco FJ, Lotz M. Autophagy is a protective mechanism in normal cartilage, and its aging-related loss is linked with cell death and osteoarthritis. Arthritis Rheum (2010) 62(3):791–801. doi: 10.1002/art.27305

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Yue J, Aobulikasimu A, Sun W, Liu S, Xie W, Sun W. Targeted regulation of foxo1 in chondrocytes prevents age-related osteoarthritis via autophagy mechanism. J Cell Mol Med (2022) 26(11):3075–82. doi: 10.1111/jcmm.17319

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Rollin R, Marco F, Jover JA, Garcia-Asenjo JA, Rodriguez L, Lopez-Duran L, et al. Early lymphocyte activation in the synovial microenvironment in patients with osteoarthritis: comparison with rheumatoid arthritis patients and healthy controls. Rheumatol Int (2008) 28(8):757–64. doi: 10.1007/s00296-008-0518-7

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Martel-Pelletier J, Barr AJ, Cicuttini FM, Conaghan PG, Cooper C, Goldring MB, et al. Osteoarthritis. Nat Rev Dis Primers (2016) 2(1):16072. doi: 10.1038/nrdp.2016.72

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Germic N, Frangez Z, Yousefi S, Simon HU. Regulation of the innate immune system by autophagy: monocytes, macrophages, dendritic cells and antigen presentation. Cell Death Differ (2019) 26(4):715–27. doi: 10.1038/s41418-019-0297-6

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Valdor R, Macian F. Autophagy and the regulation of the immune response. Pharmacol Res (2012) 66(6):475–83. doi: 10.1016/j.phrs.2012.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Huang RZ, Zheng J, Liu FL, Li QL, Huang WH, Zhang DM, et al. A novel autophagy-related marker for improved differential diagnosis of rheumatoid arthritis and osteoarthritis. Front Genet (2021) 12:743560. doi: 10.3389/fgene.2021.743560

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Fisch KM, Gamini R, Alvarez-Garcia O, Akagi R, Saito M, Muramatsu Y, et al. Identification of transcription factors responsible for dysregulated networks in human osteoarthritis cartilage by global gene expression analysis. Osteoarthritis Cartilage (2018) 26(11):1531–8. doi: 10.1016/j.joca.2018.07.012

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for rna-seq data with deseq2. Genome Biol (2014) 15(12):550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. Clusterprofiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (N Y) (2021) 2(3):100141. doi: 10.1016/j.xinn.2021.100141

CrossRef Full Text | Google Scholar

21. Walter W, Sanchez-Cabo F, Ricote M. Goplot: an R package for visually combining expression data with functional analysis. Bioinformatics (2015) 31(17):2912–4. doi: 10.1093/bioinformatics/btv300

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Lu Z, Lei D, Jiang T, Yang L, Zheng L, Zhao J. Nerve growth factor from chinese cobra venom stimulates chondrogenic differentiation of mesenchymal stem cells. Cell Death Dis (2017) 8(5):e2801. doi: 10.1038/cddis.2017.208

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Carames B, Hasegawa A, Taniguchi N, Miyaki S, Blanco FJ, Lotz M. Autophagy activation by rapamycin reduces severity of experimental osteoarthritis. Ann Rheum Dis (2012) 71(4):575–81. doi: 10.1136/annrheumdis-2011-200557

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Lopez de Figueroa P, Lotz MK, Blanco FJ, Carames B. Autophagy activation and protection from mitochondrial dysfunction in human chondrocytes. Arthritis Rheumatol (2015) 67(4):966–76. doi: 10.1002/art.39025

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Gang X, Xu H, Si L, Zhu X, Yu T, Jiang Z, et al. Treatment effect of cdkn1a on rheumatoid arthritis by mediating proliferation and invasion of fibroblast-like synoviocytes cells. Clin Exp Immunol (2018) 194(2):220–30. doi: 10.1111/cei.13161

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Huang S, Xu M, Liu L, Yang J, Wang H, Wan C, et al. Autophagy is involved in the protective effect of P21 on lps-induced cardiac dysfunction. Cell Death Dis (2020) 11(7):554. doi: 10.1038/s41419-020-02765-7

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Yang C, Xu X, Dong X, Yang B, Dong W, Luo Y, et al. Ddit3/chop promotes autophagy in chondrocytes via sirt1-akt pathway. Biochim Biophys Acta Mol Cell Res (2021) 1868(9):119074. doi: 10.1016/j.bbamcr.2021.119074

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zhang L, Xiao F, Zhang J, Wang X, Ying J, Wei G, et al. Dexmedetomidine mitigated nlrp3-mediated neuroinflammation via the ubiquitin-autophagy pathway to improve perioperative neurocognitive disorder in mice. Front Pharmacol (2021) 12:646265. doi: 10.3389/fphar.2021.646265

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Lu Y, Zhang C, Jiang S, Yuan F. Anti-dlx5 retards the progression of osteoarthritis through inhibiting chondrocyte hypertrophy and apoptosis. Evid Based Complement Alternat Med (2022) 2022:5019920. doi: 10.1155/2022/5019920

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Yuan WH, Xie QQ, Wang KP, Shen W, Feng XF, Liu Z, et al. Screening of osteoarthritis diagnostic markers based on immune-related genes and immune infiltration. Sci Rep (2021) 11(1):7032. doi: 10.1038/s41598-021-86319-7

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Deng YJ, Ren EH, Yuan WH, Zhang GZ, Wu ZL, Xie QQ. Grb10 and E2f3 as diagnostic markers of osteoarthritis and their correlation with immune infiltration. Diagnostics (Basel) (2020) 10(3):171. doi: 10.3390/diagnostics10030171

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Cai W, Li H, Zhang Y, Han G. Identification of key biomarkers and immune infiltration in the synovial tissue of osteoarthritis by bioinformatics analysis. PeerJ (2020) 8:e8390. doi: 10.7717/peerj.8390

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Han Y, Wu J, Gong Z, Zhou Y, Li H, Wang B, et al. Identification and development of a novel 5-gene diagnostic model based on immune infiltration analysis of osteoarthritis. J Transl Med (2021) 19(1):522. doi: 10.1186/s12967-021-03183-9

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Deretic V. Autophagy in inflammation, infection, and immunometabolism. Immunity (2021) 54(3):437–53. doi: 10.1016/j.immuni.2021.01.018

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Deretic V, Levine B. Autophagy balances inflammation in innate immunity. Autophagy (2018) 14(2):243–51. doi: 10.1080/15548627.2017.1402992

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Levine B, Mizushima N, Virgin HW. Autophagy in immunity and inflammation. Nature (2011) 469(7330):323–35. doi: 10.1038/nature09782

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Wang L, Fang D, Liu Y. Autophagy-related genes are potential diagnostic biomarkers for dermatomyositis. Ann Transl Med (2022) 10(4):228. doi: 10.21037/atm-22-70

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Alahdal M, Zhang H, Huang R, Sun W, Deng Z, Duan L, et al. Potential efficacy of dendritic cell immunomodulation in the treatment of osteoarthritis. Rheumatol (Oxford) (2021) 60(2):507–17. doi: 10.1093/rheumatology/keaa745

CrossRef Full Text | Google Scholar