Consensus cluster analysis of apoptosis-related genes in patients with osteoarthritis and their correlation with immune cell infiltration

Background: Osteoarthritis (OA) progression involves multiple factors, including cartilage erosion as the basic pathological mechanism of degeneration, and is closely related to chondrocyte apoptosis. To analyze the correlation between apoptosis and OA development, we selected apoptosis genes from the differentially expressed genes (DEGs) between OA and normal samples from the Gene Expression Omnibus (GEO) database, used lasso regression analysis to identify characteristic genes, and performed consensus cluster analysis to further explore the pathogenesis of this disease.

Methods: The Gene expression profile datasets of OA samples, GSE12021 and GSE55235, were downloaded from GEO. The datasets were combined and analyzed for DEGs. Apoptosis-related genes (ARGs) were collected from the GeneCards database and intersected with DEGs for apoptosis-related DEGs (ARDEGs). Least absolute shrinkage and selection operator (LASSO) regression analysis was performed to obtain characteristic genes, and a nomogram was constructed based on these genes. A consensus cluster analysis was performed to divide the patients into clusters. The immune characteristics, functional enrichment, and immune infiltration statuses of the clusters were compared. In addition, a protein–protein interaction network of mRNA drugs, mRNA-transcription factors (TFs), and mRNA-miRNAs was constructed.

Results: A total of 95 DEGs were identified, of which 47 were upregulated and 48 were downregulated, and 31 hub genes were selected as ARDEGs. LASSO regression analysis revealed nine characteristic genes: growth differentiation factor 15 (GDF15), NAMPT, TLR7, CXCL2, KLF2, REV3L, KLF9, THBD, and MTHFD2. Clusters A and B were identified, and neutrophil activation and neutrophil activation involved in the immune response were highly enriched in Cluster B, whereas protein repair and purine salvage signal pathways were enriched in Cluster A. The number of activated natural killer cells in Cluster B was significantly higher than that in Cluster A. GDF15 and KLF9 interacted with 193 and 32 TFs, respectively, and CXCL2 and REV3L interacted with 48 and 82 miRNAs, respectively.

Conclusion: ARGs could predict the occurrence of OA and may be related to different degrees of OA progression.

1 Introduction

Osteoarthritis (OA) mainly affects the articular cartilage and is the most common joint disease. It is characterized by joint pain, stiffness, hypertrophy, and limited activity. It is most common in weight-bearing joints, such as the knee, hip, cervical, and lumbar spine, as well as the proximal and distal finger joints (1). It is the most prevalent among middle-aged and older individuals. Among people aged >65 years, the incidence of the disease is as high as 50% (2). The following clinical manifestations support the diagnosis of OA: Age of onset ≥ 45 years; Persistent usage-related joint pain in one or a few joints; Early morning and Short-lived (≤ 30 minutes) stiffness (3). According to X-ray findings, OA can be classified into 5 grades as follows: 0-None; 1-Doubtful; 2-Minimal; 3-Moderate; 4-Severe (4). The differential diagnosis for OA mainly depends on the location of the affected site and the presence other systemic symptoms. The differential diagnosis includes rheumatoid arthritis, psoriatic arthritis, crystalline arthritis, hemochromatosis, infectious arthritis, and other soft-tissue abnormalities. This disease is a progressive disease that eventually causes degeneration, fibrosis, rupture, defects, and damage to the entire articular surface of the joint, considerably affecting the lives of patients. If not treated promptly, there may be a risk of disability for patients (5). Owing to the fact that the etiology and pathogenesis of OA are not yet fully understood, symptomatic treatment can only be used to delay the progression of the disease, alleviate the symptoms of patients, correct deformities through surgery in the late stage, and improve the function of affected limbs (6). Therefore, OA pathogenesis requires further investigation.

Currently, most researchers believe that OA is the result of a combination of mechanical and biological factors leading to an imbalance in the degradation and anabolism of articular chondrocytes, extracellular matrix (ECM), and subchondral bone (7). Apoptosis is a process of programmed cell death controlled by genes. The human body maintains the homeostasis of various system tissues by clearing the dead cells and metabolites. However, if the apoptosis process is disrupted, it is directly or indirectly related to the occurrence of many diseases, such as tumors, autoimmune diseases, local injuries, etc. (8). Apoptosis in the cartilage tissue related to matrix degradation and calcification has been detected in OA cartilage, indicating that apoptosis plays a role in OA pathogenesis (9).

To comprehensively understand the role of apoptosis and immune mechanisms in the occurrence and development of OA, we analyzed transcriptome data from the GEO database, selected characteristic apoptosis-related genes (ARGs) that could predict OA occurrence, and attempted to cluster patients with OA based on these genes. We analyzed the immune characteristics, biological pathways, and immune infiltration between different clusters and searched for drugs and molecules that inhibit OA occurrence.

2 Materials and methods

2.1 Data downloading and initial preparation

OA-related chip sequencing data and corresponding clinical sample information were downloaded from the GSE12021 (10) and GSE55235 (11) datasets in the GEO database. The sample source was Homo sapiens and the two dataset sequencing platforms were GPL96 and GPL97. The data set GSE12021 contains 33 samples in total, including 13 normal samples and 20 samples from patients with OA; the data set GSE55235 contains 20 samples in total, including 10 normal samples and 10 samples from patients with OA. We integrated the above sets of data expression into one data set, which was named as the combined data or combined gene expression. Then, the R package “sva” (12) was used to correct the batch effect between different data sets and log2 (X+1) standardization was performed. Standardized data were used as the basis for subsequent analyses.

2.2 Apoptosis-related differentially expressed genes

A total of 3,576 ARGs were collected from the GeneCards database (13) and other documents. The differential gene analysis between OA samples and normal samples was conducted using the R package “limma” (14). Significant DEGs were screened using the absolute value of the log2Fold change (log2FC) > 1 and P< 0.05. The upregulated genes were defined as log2FC > 1 and P< 0.05, and the downregulated DEGs were defined as log2FC< - 1 and P< 0.05. The DEG expression results are displayed on a volcano map and a heatmap. The ARGs and DEGs are intersected as hub genes called ARDEGs.

2.3 Least absolute shrinkage and selection operator regression analysis and risk model construction

To accurately screen for biomarkers related to OA, we conducted dimensionality reduction screening using the LASSO model with 1,000 iterations (15). The objective function of the LASSO regression model was as follows:

$$min{\int }^{ }({\alpha }_0,\alpha |X_i,Y_i+ \lambda {\left|\left|\alpha \right|\right|}_1)$$

where λ represents the penalty coefficient, which can be selected through 10-fold cross-validation for the chosen λ; ||α||₁ is defined as the sum of the absolute values of each vector element. LASSO regression was implemented through the R package “glmnet” (16). A risk-scoring formula was established by weighing each normalized gene expression value with the penalty coefficient of the characteristic gene.

$$\text{riskScore }=\sum _{\text{i}}\text{Coefficient }({\text{gene}}_{\text{i}})*\text{mRNA Expression }({\text{ gene}}_{\text{i}})$$

Subsequently, a nomogram model was constructed based on selected candidate ARDEGs to predict OA.

2.4 Consensus clustering method and identification of apoptotic subtypes

Based on the expression data of the characteristic genes in the Combine data chipset with all the samples, apoptotic clusters were identified using the R package “ConsensusClusterPlus” (17). The R package “limma” was used again to analyze the differential genes of different apoptotic subtypes in the combined data set, with |log2FC| > 1 and P< 0.05 as the differential gene-screening criteria. DEGs with |log2FC| > 4 and P< 0.05 were considered upregulated, whereas those with |log2FC|< -3 and P< 0.05 were considered downregulated. The results are displayed using a volcano map.

2.5 Gene function enrichment analysis by GO and single-sample gene set enrichment analysis

GO enrichment analysis (18) is a common method for conducting large-scale functional enrichment studies on biological processes (BPs), molecular functions (MFs), and cellular components (CCs). It was conducted for all the ARDEGs through the “clusterProfiler” package (19). To analyze the differences in BPs between subgroups based on gene expression data, we used ssGSEA (20), which is a computational method that analyzes whether a specific gene set shows statistical differences between two biological states. It is commonly used to estimate changes in pathways and BPs in expression dataset samples, with a P value< 0.05 considered to be statistically significant.

2.6 Correlation analysis and chromosome location between genes

To explore the correlation among ARDEGs, spearman analysis was performed through the R package “cowplot” (21), and heat maps, scatter maps, and correlation curves were drawn. P values< 0.05 imply that the genes have a strong correlation. The R package “RCircos” (22) was used to draw the location map of the key genes in the chromosome with the location information of genes downloaded from the ENSEMBL database (23).

2.7 CIBERSORT for the infiltration state of immune cells

CIBERSORT (https://cibersort.stanford.edu/) (24) is an R/webpage version tool that deconvolutes the expression matrix of human immune cell subtypes based on the principle of linear support vector regression. The expression of 22 types of known immune cells was calculated and evaluated for infiltration status using the CIBERSORT algorithm, and the percentage difference of immune cells between the OA samples and the normal samples was tested using the Wilcoxon test, with a P value< 0.05 considered to be statistically significant.

2.8 Drug–gene interaction analysis

We searched the DGIdb database (Version 3.0.2, https://www.dgidb.org) (25) with ARDEG names as the input to predict any potential drugs or molecular compounds interacting with them. The results were visualized by a drug–gene interaction network using the Cytoscape software.

2.9 Construction of mRNA–miRNA and mRNA–transcription factor networks

MiRNAs (26) are a class of non-coding single-stranded RNA molecules encoded by endogenous genes, with a length of approximately 19–25 nucleotides, that play an essential regulatory role in the evolution of biological development. MiRNAs regulate the expression of target genes by participating in the post-transcriptional regulation of genes and play a crucial regulatory role in tumor occurrence and development, biological development, organ formation, epigenetic regulation, and viral defense. Typically, miRNAs have very complex regulatory networks and often one miRNA can regulate multiple target genes, whereas the same target gene can also regulate many miRNAs. To analyze the relationship between hub genes and miRNAs at the post-transcriptional stage, we obtained hub gene miRNAs from the Network Analyst database (27) and constructed an mRNA–miRNA regulatory network. TFs control gene expression by interacting with target genes during the post-transcriptional stage. To analyze the regulatory effect of TFs on hub genes, the targeted relationship between transcription factors and hub genes was retrieved from the Network Analyst database, and an interaction network between hub genes and TFs was constructed.

2.10 Statistical methods

All data were processed and analyzed using the R software (version 4.1.1). To compare two continuous variables that were distributed normally, an independent Student’s t-test was conducted. The Mann–Whitney U test was performed to determine the difference between non-normally distributed variables. Pearson’s correlation analysis was used to calculate the correlation coefficient (r) between the different genes, the Pearson’s correlation coefficient between two variables is defined as the quotient of the covariance and standard deviation between two variables:

$$\rho X,Y= \frac{cov(X,Y)}{{\sigma }_X{\sigma }_Y}= \frac{E[(X-{\mu }_X)(Y-{\mu }_Y)]}{{\sigma }_X{\sigma }_Y}$$

The above equation defines the overall correlation coefficient, and the Greek lowercase letters are commonly used as representative symbols. To estimate the covariance and standard deviation of the sample, the Pearson correlation coefficient (r). P values of all statistical processes were bilateral, with P values< 0.05 considered significant.

3 Results

3.1 Workflow

The workflow is depicted in Figure 1.

FIGURE 1

Figure 1 The flow chart.

3.2 Data set combination, OA-related DEGs, and their GO functional enrichment

In this study, the two data sets GSE12021 and GSE55235 were combined, and the batch effect between the data sets was removed (Figure 2). Ninety-five DEGs were then obtained from the OA (n = 30) and control samples (n = 23), including 47 upregulated and 48 downregulated genes (Figures 3A, B). GO functional annotation was performed on the DEGs to determine their biological functions (Figures 3C–E). The results showed that these DEGs were mainly enriched in BPs, such as ECM organization, extracellular structure organization, and CCs, such as collagen-containing extracellular matrix, immunoglobulin complex/circulating, as well as MFs such as immunoglobulin receptor binding, antigen binding (Figures 3C–E).

FIGURE 2

Figure 2 Gene expression distribution and batch effect correction of all samples. (A) The difference between the two data sets before removing the batch effect. (B) The difference between the two data sets after removing the batch effect.

FIGURE 3

Figure 3 The ARDEGs. (A) The volcano map of DEGs related to osteoarthritis (OA). The x-axis represents log2FoldChange, and the y-axis represents -log10 (P-value). The red node represents the upregulated differentially expressed genes (DEGs), blue node represents the downregulated DEGs, and the gray node represents the non-significant genes. (B) Heat map of the expression level of DEGs related to OA; green represents disease samples, brown represents normal control samples, red represents high expression, and blue represents low expression. (C–E) Functional enrichment analysis of DEGs. (C) The abscissa represents the enriched gene book, the ordinate represents the biological process, and the color indicates the significance of the enrichment results. (D) Ring diagram of GO enrichment results. The innermost ring in the figure is the gene cluster tree. The middle ring is the logFC value corresponding to these genes. The closer the red logFC value, the higher the upregulation level. The outer ring is the GO terms enriched by these genes. Each GO term uses a color to set itself apart. (E) The result of the enrichment analysis of the KEGG pathway shows that the node color represents the gene expression level, and the quadrilateral color represents the KEGG pathway Z-core. ARDEGs, apoptosis-related differentially expressed genes; DEGs, differentially expressed genes.

3.3 Risk prediction model construction

A total of 31 hub genes were selected from the apoptotic genes and DEGs (Figure 4A), and hub genes with FDR< 0.05 were selected as candidate genes. The hub genes were put into logistic Lasso analysis. The analysis results revealed nine characteristic genes, namely growth differentiation factor 15 (GDF15), NAMPT, TLR7, CXCL2, KLF2, REV3L, KLF9, THBD, and MTHFD2, which occurred 413 times in 1,000 cycles (Figure 4B; Table 1). The gene correlations of the nine characteristic genes were calculated. The results showed a strong correlation between CXCL2 and KLF2 (r = 0.8, P< 0.05, r for the correlation coefficient), also CXCL2 and THBD (r = 0.74, P< 0.05) (Figure 4C). The location of 18 genes on the human chromosome through the R package RCircos is shown (Figure 4D), and the results showed that these genes were frequently found on chr2, chr4, chr6, chr7, chr9, chr19, chr20, and chrX.

TABLE 1

Table 1 The gene lists.

FIGURE 4

Figure 4 Construction and correlation analysis of the diagnosis model of osteoarthritis. (A, B) Through lasso algorithm, characteristic genes are selected from apoptosis-related genes. The horizontal axis is the combination of genes selected through 1000 Lasso analyses, and the vertical axis is the number of occurrences. (C) For the correlation analysis of nine characteristic genes in all samples, *P< 0.05, **P< 0.01, ***P< 0.001, and the number represents the correlation level. (D) Location map of nine characteristic genes on the chromosome.

Based on the patient’s predicted risk score for OA, the results showed that the score had a profound effect (Figure 5A). The calibration curve indicated that the nomogram model was accurate (Figure 5B). Decision curve analysis (DCA) was used to evaluate the potential clinical impact of the nomogram model and its association with OA. The dashed line in the DCA curve remained above the grey and black lines from 0 to 1, indicating that the decision based on the nomogram model may be beneficial for patients with OA (Figure 5C). Simultaneously, the ROC curves of the nine characteristic genes were analyzed to predict OA independently, and the results showed that the nine characteristic genes had an efficient predictive ability (Figure 5D).

FIGURE 5

Figure 5 Nomogram. (A) Nomograms predicting the risk score for the diagnosis of patients with osteoarthritis (OA). (B) The model evaluation curve; gray represents random diagnosis, while green and pink represent model diagnosis. (C) DCA curve. The x-axis represents the risk threshold, y-axis represents the net profit rate, black solid line represents 0 net profit rate, all (gray solid line) represents that all samples are subject to intervention, and the affected site modeling (black dotted line) represents the model curve. (D) ROC curve of nine characteristic genes in the diagnosis of OA.

3.4 Apoptosis mode identification

Using the “ConsensusClusterPlus” package in the R software, based on nine characteristic genes, two cell apoptosis subtypes (ClusterA and ClusterB) were identified using the consistency clustering method (Figures 6A, B). The volcanic map of DEGs of Cluster A and Cluster B was shown in the Figure 6C. Except for TLR7, the nine characteristic genes were expressed at high levels in the normal samples. ClusterA contained 13 samples and ClusterB contained 17 samples. The expression levels of the characteristic genes CXCL2, GDF15, KLF2, MTHFD2, NAMPT, and REV3L were significantly higher in ClusterA than those in ClusterB, whereas the expression level of KLF9 was higher in ClusterB than that in ClusterA (Figures 6D, E).

FIGURE 6

Figure 6 Consistency cluster score of patients with osteoarthritis (OA) and functional analysis between the two apoptosis modes. (A) Consistency clustering result diagram. (B) CDF, this figure shows the cumulative distribution function when K takes different values. It is used to judge the value of K, and CDF reaches the approximate maximum value. During this time, the cluster analysis results are the most reliable, usually taking the value of K with a small decline slope of CDF. (C) Cluster A and Cluster B differentially expressed genes (DEGs) have a volcanic map with log2FoldChange as the abscissa and -log10 (P-value) as the ordinate. The red node represents the upregulated DEGs, green node represents the downregulated DEGs, and gray node represents the cause of non-significant DEGs. (D) Histograms of the expression levels of nine characteristic genes in the disease group and the control group, with blue indicating the control group and yellow indicating the disease group. (E) Histogram of the expression level of nine characteristic genes in ClusterA and ClusterB of the OA subgroups, blue represents ClusterA, and yellow represents ClusterB. (F) GSEA analysis results between ClusterA and ClusterB. (G) SSGSEA analysis of Cluster A and Cluster B, KEGG biological function enrichment analysis. CDF, cumulative distribution function; SSGSEA, single sample gene set enrichment analysis. * represents P < 0.1, ** represents P < 0.01, *** represents P < 0.001, ns represents non-significant.

3.5 Analysis of functional differences between two different cell apoptosis modes

To analyze the differences between the two different cell apoptosis modes, DEGs were obtained in the two cell apoptosis modes, ClusterA and ClusterB (Figure 6F). We then analyzed the effects of DEGs between the two cell apoptosis modes on the biological functions of patients. GSEA revealed that neutrophil activation and neutrophil activation involved in immune response were highly enriched in ClusterB (Figure 6F). Using ssGSEA, we explored the highly enriched biological signaling pathways in the two cell apoptosis subgroups. The results showed that the protein repair and purine salvage signaling pathways were highly enriched in the ClusterA model (Figure 6G).

3.6 Differences in immune characteristics between the two cell apoptosis modes

The results of the CIBERSORT analysis showed that the level of activated natural killer (NK) cells in ClusterB was significantly higher than that in ClusterA. The M2 macrophage content in patients in ClusterA was significantly higher than that in ClusterB (Figures 7A, B).

FIGURE 7

Figure 7 Immune characteristics between the two apoptosis modes - CIBERSORT. (A) Accumulation diagram of immune cell content between Cluster A and Cluster B shows different immune cells in different colors, and the x-axis represents the cell fraction and the y-axis represents the sample name. (B) Histogram of the content of immune cells in patients of Cluster A and Cluster B; blue represents the Cluster A sample, and yellow represents the Cluster B sample. (C) Correlation between nine characteristic genes and immune cells. The x-axis represents the amount of gene expression, and the y-axis represents score of immune cell infiltration. r > 0 represents a positive correlation and r< 0 represents a negative correlation. * represents P < 0.1, ** represents P < 0.01, *** represents P < 0.001, ns represents non-significant.

3.7 Analysis of correlation between key genes and immune cell infiltration

The results showed that the expression of CXCL2, GDF15, KLF2, and MTHFD2 genes was positively correlated with activated NK cells (r = 0.37, P = 0.04; r = 0.68, P< 0.01; r = 0.65, P< 0.01; R = 0.56, P< 0.01, respectively), whereas the gene expression of KLF9 was negatively correlated with activated NK cells (r = -0.6, P< 0.01) (Figure 7C). The expression of CXCL2, GDF15, KLF2, and MTHFD2 genes was negatively correlated with M2 macrophages (r = -0.45, P = 0.01; r = -0.8, P< 0.01; r = -0.63, P< 0.01; r = -0.53, P< 0.01), whereas the gene expression of KLF9 was positively correlated with M2 macrophages (r = 0.6, P< 0.01) (Figure 7C).

The ssGSEA results showed that the contents of several immune cells in ClusterA were higher than those in ClusterB, such as the activated CD8 T cell, the central memory CD4 T cell, the central memory CD8 T cell, and the Type 1 T helper cell (P< 0.05). Nevertheless, the contents of some immune cells in ClusterA were lower than those in ClusterB, such as regulatory T cell, the natural killer cell, the gamma delta T cell, Type 2 T helper cell, the neutrophil, and nature killer T cell (P< 0.05, Figure 8A). Correlation analysis showed that the nine characteristic genes were closely related to immune cell content (Figure 8B).

FIGURE 8

Figure 8 Immune characteristics between the two cell apoptosis modes - ssGSEA. (A) Histogram of the content of immune cells in patients of Cluster A and Cluster B; yellow represents the Cluster B sample and blue represents the Cluster A sample. (B) Correlation analysis between the nine characteristic genes and the content of immune cells shows that the x-axis represents the immune cells, y-axis represents the characteristic genes, node color represents the correlation size, and node size represents the significance level. ssGEA, single-sample gene set enrichment analysis. * represents P < 0.1, ** represents P < 0.01, *** represents P < 0.001, ns represents non-significant.

As shown in Figure 9, the correlations between the immune cells in patients with OA was further calculated, the results showed that NK cells showed significantly positively correlated with gamma delta T cells, MDSC, and eosinophils (P< 0.05), whereas negatively correlated with immature B cells, macrophages, and activated dendritic cells (P< 0.05).3.8 Drug–gene interaction analysis.

FIGURE 9

Figure 9 Correlation between immune cells in patients with osteoarthritis (OA). The correlation of immune cells in patients with OA is negative in red and positive in blue.

The DGIdb database was used to retrieve small-molecule compounds and drugs that regulated key genes. As shown in the drug–gene interaction network, 29 drugs or molecular compounds, including mepyrazole, conafinil, and abacavir, were associated with the key genes NFKBIA, ARAF, and ADH1B of OA ClusterA (Figure 10A). A total of 38 drugs or molecular compounds, including PICTILISIB, Linifanib, and CHEMBL225519, were related to the key genes STK11, PLK3, and PADI1 of OA ClusterB (Figure 10B). The relevance of these drugs or compounds to key genes could imply that they may have varying degrees of effect in regulating these genes.

FIGURE 10

Figure 10 Drug sensitivity analysis. (A) Cell apoptosis of osteoarthritis is related to small molecular compounds or drug sensitivity analysis of Cluster A. (B) Cell apoptosis of osteoarthritis is related to small molecular compounds or drug sensitivity analysis of Cluster B. Red represents the core differential genes of each subtype, gray represents small molecular compounds or drugs, and the connection represents the connection between drugs and genes.

3.9 PPI network

We constructed mRNA–miRNA and mRNA–TF networks of the characteristic genes. The mRNA-TF network included 1,324 interactions, 8 mRNA, and 368 TFs. GDF15 interacted with 193 TFs, while KLF9 interacted with 232 TFs (Figure 11A). The mRNA–miRNA network of the nine characteristic genes was constructed, including 4,213 interaction relationships, 9 mRNAs, and 1,236 miRNAs, with CXCL2 interacting with 48 miRNAs and REV3L interacting with 82 miRNAs (Figure 11B).

FIGURE 11

Figure 11 mRNA–miRNA and mRNA–TF network associated with characteristic genes. (A) The mRNA–TF network related to characteristic genes, with gray nodes representing TF and red nodes representing characteristic genes. (B) The mRNA-miRNA network associated with the characteristic gene. The gray node represents the miRNA and the blue node represents the characteristic gene. TF, transcription factor.

4 Discussion

OA is the most common chronic degenerative joint disease encountered in orthopedic clinics, among which the most commonly affected joints are the hip and knee joints. The incidence rate of OA increases annually with increasing age, especially in middle-aged and older populations (28). In addition to relieving symptoms, there are currently no reports confirming that existing treatment measures can prevent or reverse OA progression. Therefore, the pathogenesis of atherosclerosis has become a popular research topic. OA is the result of multiple factors, such as cartilage nutrition, metabolic abnormalities, stress imbalance, abnormal degradation of the cartilage matrix by enzymes, cumulative minor trauma, obesity, and increased joint weight bearing (29). The occurrence of OA is a continuous process, and its basic pathological mechanism involves cartilage degeneration and erosion. As the disease continues to develop, cystic changes in the subchondral bone occur gradually, and bone spurs begin to appear at the edge of the joint. In this continuous pathological change, cartilage degeneration is a key link and the initiating factor in the occurrence and development of OA (30). In addition to relieving the symptoms of OA, there are currently no reports confirming whether the existing treatment measures can prevent or reverse its progression. Therefore, studying OA pathogenesis and progression is essential for understanding the disease and proposing effective treatment targets.

Some researchers believe that OA development is closely related to chondrocyte apoptosis. Hwang et al. (31) suggested that significant apoptosis occurs in advanced OA, resulting in a decrease in chondrocytes, often accompanied by lacunar emptying, highlighting that chondrocyte apoptosis is a characteristic of OA progression. However, the role of apoptosis in the progressive degeneration of cartilage remains unclear. Identifying key genes and pathways related to cell apoptosis can help people understand the process of OA occurrence and development and carry out targeted molecular research to design targeted therapeutic drugs.

Therefore, to analyze the correlation between apoptosis and OA development, we selected apoptosis genes from the DEGs between OA and normal samples and performed LASSO regression analysis to identify nine characteristic genes: GDF15, NAMPT, TLR7, CXCL2, KLF2, REV3L, KLF9, THBD, and MTHFD2. We also constructed a risk-prediction nomogram. Testing the prediction accuracy of the nomogram showed that nine apoptotic genes accurately predicted OA occurrence. This indicated that apoptosis and apoptotic genes play key roles in OA occurrence and development.

GDF15, a member of the tumor growth factor superfamily, has recently been identified as a possible biomarker of aging and is associated with various clinical conditions, including coronary artery disease, diabetes, and various cancer types. Wen et al. (32) studied the role of the GDF15/MAPK14 axis in the aging of chondrocytes in OA and found that GDF15 is a driving factor for chondrocyte aging and apoptosis and could promote the progression of OA by inducing angiogenesis. TLR7 is involved in immune responses in many inflammatory diseases. Liu et al. (33) believe that silencing TLR7 could block the p21-mediated JAK2/STAT3 pathway and prevent lipopolysaccharide-induced apoptosis and chondrocyte injury. The results of the present study are consistent with those of these studies.

We conducted a differential analysis using OA samples and control samples to obtain 95 DEGs, which were mainly enriched in ECM organization, extracellular structure organization, collagen-containing ECM, immunoglobulin complex, circulation, immunoglobulin circulating, immunoglobulin-receiver-binding, and antigen binding. Some researchers have suggested that with the progression of the disease, the immune barrier of the ECM disappears and specific surface antigens of chondrocytes cause autoimmune reactions, opening a molecular damage mode for ECM products. The innate immune and inflammatory cycle mechanisms in OA lead to sustained joint damage (34). The results of our genetic enrichment analysis are consistent with these findings.

Previous studies have found rare mutations in certain genes that can lead to hereditary OA, even results in severe OA (35–37). At least 100 risk loci associated with OA have also been identified by genome wide-association studies (38, 39). In this study, we plotted the genetic locus on the chromosome of characteristic genes. These results provide clues for further studies on the genetic underpinnings and biological mechanisms of OA pathogenesis. First, neighboring genes in the same chromosomal regions as the potentially OA-related genes may participate in similar pathways and functions. Second, further studies will focus on exploring the relationships between single nucleotide polymorphisms in these genes and OA risk. Furthermore, the overlap between gene locations and regulatory elements suggests these genes could be co-regulated in certain biological processes important to OA. In our study, differential expression gene analysis revealed interesting patterns, such as KLF2 and GDF15 genes located on the same chromosome. The observation of differential genes clustered on specific chromosomes has prompted us to delve deeper into the potential implications of this phenomenon. The co-localization of differentially expressed genes on the same chromosome may reflect a coordinated regulatory mechanism that regulate their expression levels. This co-regulation may be attributed to various factors, such as shared transcriptional regulatory elements, epigenetic modifications, and even physical chromatin interactions. By enhancing our understanding of OA genetics and biology, these findings could ultimately promote drug development.

To subdivide the different pathways and characteristics of cell apoptosis, we attempted to use nine genes to perform a consensus cluster analysis on 30 OA samples and finally obtained two clusters, namely Cluster A and Cluster B, including 13 samples for ClusterA and 17 samples for Cluster B. The expression of nine characteristic genes in both clusters was significantly increased or decreased. We compared the biologically related functions and pathways of the two clusters and found that cluster A was highly enriched in protein repair and purine salvage signaling pathways, whereas Cluster B was highly enriched in neutrophil activation, neutrophil activation involved in immune response, and other neutrophil-related pathways in the cluster B subgroup. These results suggest that cluster A was highly enriched in cell repair and rescue, and cluster B was involved in inflammation-related immune responses. Immune cell infiltration analysis showed that the expression levels of T cells and M2 macrophages was higher in cluster A, whereas the expression levels of activated NK cells was significantly higher in cluster B. This reflects the different immune microenvironments of the two types of samples. We speculated that this difference in immune cell infiltration may be related to the different OA progression in the two types of samples. Chondrocyte aging and apoptosis play a promoting role during OA occurrence and development, (Cluster B). At this stage, the patient’s immune response is dominated by adaptive immunity, manifested by the disappearance of the ECM immune barrier and the release of specific surface antigens from chondrocytes, which stimulate the corresponding T cells in the lymphocyte pool of the immune system, triggering the process of cell-mediated immunity. In this process, T cells release lymphatic factors to clear antigens, causing damage to the cartilage cells and initiating BPs, such as protein repair and purine rescue. However, the high level of NK cells is related to the clearance of inflammatory factors, and the high concentration of neutrophil-related pathways in cluster A suggests that the patients have progressed to the stage of inflammatory circulation, and its immune response is dominated by innate immunity, with an increase in the level of NK cells as a significant feature.

We also observed that the expression levels of the CXCL2, GDF15, KLF2, and MTHFD2 genes were positively correlated with the number of activated NK cells and negatively correlated with the number of M2 macrophages; however, the reverse was true for KLF9. Therefore, we speculate that these five genes are involved in regulating the immune response. Evidence for this can be found in various literature. CXCL2 is a member of the chemokine superfamily that encodes secreted proteins involved in immune regulation and inflammatory processes. This chemokine is a member of the CXC subfamily, expressed at inflammatory sites, and may inhibit the proliferation of hematopoietic progenitor cells (40). KLF2 regulates innate immune responses during skeletal muscle injury and regeneration (41). Wang et al. (42) confirmed that GDF15 induces immunosuppression through regulatory T cells.

In colon adenocarcinoma, patients with high GDF15 had favorable overall survival, but this survival advantage was reversed when they also had decreased NK cells (43). Another study suggested that decidual NK cells contribute to embryo growth by GDF15 secretion (44). Besides, GDF15 overexpression was associated with increased CD8 T cell numbers and proportion of activated CD8+CD11c+ T cells. Furthermore, depletion of CD8 T cells in tumor-bearing mice eliminated the protective effect of GDF15 against tumor growth (45). In accordance with these observations, our results suggested that both GDF15 and NK cells and activated CD8 T cells were significantly increased in Cluster A, and they were positively correlated. Further investigations with some wet lab evidences are needed to validate our data and find more meaningful results.

According to our results, there were many changes and correlations among immune cell subpopulations, which indicated that various immune cells contribute to pathogenesis of OA through complex mechanisms. Macrophages, B cells, T cells, and NK cells produced proinflammatory cytokines like interleukin-1β (IL-1β) (46). which induces extracellular matrix degradation (47) and chondrocyte apoptosis by promoting chondrocyte hypertrophy and dedifferentiation (31). It has been reported that Th17 cells drive OA progression by stimulating osteoclast progenitor recruitment through increased chemokine production from bone marrow mesenchymal stromal cells (48). In the patients with OA, the number of CD4+ T cell is elevated in the subsynovial layer compared to healthy group (49). CD4+ T cell secret interferon-gamma (IFN-γ) which can promote mesenchymal stem cell differentiation, while transforming growth factor-β (TGF-β) from these cells is negatively related to osteoblast differentiation (50). Distinct CD4+ T helper cell subsets also influence joint inflammation and remodel. Th1 cells promote immune response via producing IFN-γ- and IL-2-related cytokines (51). IL-4 released by Th2 cells exerts protective effects on cartilage (49). Additionally, memory T cells impair immune function (52). Targeting these various immune cell types and their signaling molecules may be helpful to clarify the pathogenesis of OA and explore new therapeutic strategies for managing OA.

Treg cells have been found to participate in changes in subchondral bone remodeling, an important process in the pathogenesis of OA. Treg cells secret various of cytokines and further promote osteoclast maturation (53). However, Treg cells also secrete IL-17F, IL-17A and BMP-2 that strongly promote osteoblastic differentiation (54). Additionally, Treg cells can inhibit bone resorption by promoting osteoclast apoptosis (55). While Treg cells have complex effects on bone remodeling through multiple mechanisms, the exactly biological mechanisms require further exploration.

T helper (Th) cells interact with other immune cells, such as B lymphocytes, through cytokine signaling and participate in immune regulation (56). B lymphocytes can activate T cells by presenting antigens like granulocyte colony-stimulating factor (G-CSF) (57), which promotes osteoclast progenitor proliferation (58). Inflammatory mediators secreted by neutrophils and osteoclasts differentially impact mesenchymal stem cells (MSCs) (59). MSCs can inhibit NK cell cytotoxicity and proliferation, prevent autoreactive antibody production, suppress Th1 cell activation, and stimulate Treg cell generation (60). Additionally, macrophages communicate with osteocytes through paracrine signaling and direct cell-cell contact (61). Taken together, interactions among various immune cell types including T cells, B cells, neutrophils, osteoclasts, MSCs, and macrophages influence bone remodeling and the progression of OA.

In drug sensitivity analysis, we compared the differential genes of two modes of cell apoptosis (Cluster A and Cluster B) and attempted to find targeted therapeutic drugs. We found 29 molecular compounds targeting key genes of Cluster A and 38 molecular compounds targeting Cluster B. Although these drugs were rarely used in the clinical treatment of osteoarthritis, they may become potential treatment agents, which remain to be discovered in the future studies.

This study has some limitations. First, although there are significant differences in the immune characteristics and functional enrichment between Cluster A and Cluster B, we have not yet found an exact grade of OA to correspond to the clusters; second, our data were from GEO. The total sample size for OA was only 30 cases, which is relatively small, and a data bias may have occurred during the research process. Third, we did not directly identify specific drugs or molecules that regulate the nine characteristic genes that affect OA progression, which should be further explored in future studies.

Apoptotic genes play a key role in the development and occurrence of OA. The ARGs GDF15, NAMPT, TLR7, CXCL2, KLF2, REV3L, KLF9, THBD, and MTHFD2 could independently and accurately predict OA occurrence. Patients with OA can be divided into two clusters with significant differences in immune characteristics, functional enrichment, and immune infiltration, which may be related to different degrees of OA progression. However, more clinical evidence is needed to confirm this. CXCL2, GDF15, KLF2, MTHFD2, and KLF9 are strongly correlated with immune infiltration in patients with OA and could become novel therapeutic targets that affect OA progression.

Data availability statement

The original contributions presented in the study are publicly available. This data can be found here: https://www.jianguoyun.com/p/DbKXe40Q1ujHCxig64IFIAA.

Author contributions

EY and JY conceived and designed the research; EY, MZ, GX and XL performed the research, analyzed the data, and prepared figures; JY wrote the paper. Al authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

Funding from the National Natural Science Foundation of China. Project approval number: 82072472.

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.1202758/full#supplementary-material

References

1. Abramoff B, Caldera FE. Osteoarthritis: pathology, diagnosis, and treatment options. Med Clin North Am (2020) 104(2):293–311. doi: 10.1016/j.mcna.2019.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Sun X, Zhen X, Hu X, Li Y, Gu S, Gu Y, et al. Osteoarthritis in the middle-aged and elderly in China: prevalence and influencing factors. Int J Environ Res Public Health (2019) 16(23):4701. doi: 10.3390/ijerph16234701

PubMed Abstract | CrossRef Full Text | Google Scholar

3. National Collaborating Centre for Chronic Conditions (UK). Osteoarthritis: National Clinical Guideline for Care and Management in Adults. London: National Institute for Health and Clinical Excellence: Guidance (2008).

Google Scholar

4. Kellgren JH, Lawrence JS. Radiological assessment of osteo-arthrosis. Ann Rheum Dis (1957) 16(4):494–502. doi: 10.1136/ard.16.4.494

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Mahmoudian A, Lohmander LS, Mobasheri A, Englund M, Luyten FP. Early-stage symptomatic osteoarthritis of the knee - time for action. Nat Rev Rheumatol (2021) 17(10):621–32. doi: 10.1038/s41584-021-00673-4

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Xia B, Di C, Zhang J, Hu S, Jin H, Tong P. Osteoarthritis pathogenesis: A review of molecular mechanisms. Calcif Tissue Int (2014) 95(6):495–505. doi: 10.1007/s00223-014-9917-9

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Prieto-Alhambra D, Judge A, Javaid MK, Cooper C, Diez-Perez A, Arden NK. Incidence and risk factors for clinically diagnosed knee, hip and hand osteoarthritis: influences of age, gender and osteoarthritis affecting other joints. Ann Rheum Dis (2014) 73(9):1659–64. doi: 10.1136/annrheumdis-2013-203355

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Xu X, Lai Y, Hua ZC. Apoptosis and apoptotic body: disease message and therapeutic target potentials. Biosci Rep (2019) 39(1):BSR20180992. doi: 10.1042/BSR20180992

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Woodell-May JE, Sommerfeld SD. Role of inflammation and the immune system in the progression of osteoarthritis. J Orthop Res (2020) 38(2):253–7. doi: 10.1002/jor.24457

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Huber R, Hummert C, Gausmann U, Pohlers D, Koczan D, Guthke R, et al. Identification of intra-group, inter-individual, and gene-specific variances in mrna expression profiles in the rheumatoid arthritis synovial membrane. Arthritis Res Ther (2008) 10(4):R98. doi: 10.1186/ar2485

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Woetzel D, Huber R, Kupfer P, Pohlers D, Pfaff M, Driesch D, et al. Identification of rheumatoid arthritis and osteoarthritis patients by transcriptome-based rule set generation. Arthritis Res Ther (2014) 16(2):R84. doi: 10.1186/ar4526

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics (2012) 28(6):882–3. doi: 10.1093/bioinformatics/bts034

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Safran M, Dalah I, Alexander J, Rosen N, Iny Stein T, Shmoish M, et al. Genecards version 3: the human gene integrator. Database (Oxford) (2010) 2010:baq020. doi: 10.1093/database/baq020

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for rna-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Wu L, Zhou B, Liu D, Wang L, Zhang X, Xu L, et al. Lasso regression-based diagnosis of acute st-segment elevation myocardial infarction (Stemi) on electrocardiogram (Ecg). J Clin Med (2022) 11(18):5408. doi: 10.3390/jcm11185408

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Ueno D, Kawabe H, Yamasaki S, Demura T, Kato K. Feature selection for rna cleavage efficiency at specific sites using the lasso regression model in arabidopsis thaliana. BMC Bioinf (2021) 22(1):380. doi: 10.1186/s12859-021-04291-5

CrossRef Full Text | Google Scholar

17. Wilkerson MD, Hayes DN. Consensusclusterplus: A class discovery tool with confidence assessments and item tracking. Bioinformatics (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U.S.A. (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Wilke CO. cowplot: Streamlined Plot Theme and Plot Annotations for 'ggplot2' (2016). Available at: https://CRAN.R-project.org/package=cowplot.

Google Scholar

20. Hanzelmann S, Castelo R, Guinney J. Gsva: Gene set variation analysis for microarray and rna-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

21. Wilke CO, Sawyer SL. At the mercy of viruses. Elife (2016) 5:e16758 doi: 10.7554/eLife.16758

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Zhang H, Meltzer P, Davis S. Rcircos: an R package for circos 2d track plots. BMC Bioinf (2013) 14:244. doi: 10.1186/1471-2105-14-244

CrossRef Full Text | Google Scholar

23. Howe KL, Achuthan P, Allen J, Allen J, Alvarez-Jarreta J, Amode MR, et al. Ensembl 2021. Nucleic Acids Res (2021) 49(D1):D884–D91. doi: 10.1093/nar/gkaa942

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol (2019) 37(7):773–82. doi: 10.1038/s41587-019-0114-2

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Cotto KC, Wagner AH, Feng YY, Kiwala S, Coffman AC, Spies G, et al. Dgidb 3.0: A redesign and expansion of the drug-gene interaction database. Nucleic Acids Res (2018) 46(D1):D1068–D73. doi: 10.1093/nar/gkx1143

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Lu TX, Rothenberg ME. Microrna. J Allergy Clin Immunol (2018) 141(4):1202–7. doi: 10.1016/j.jaci.2017.08.034

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Zhou G, Soufan O, Ewald J, Hancock REW, Basu N, Xia J. Networkanalyst 3.0: A visual analytics platform for comprehensive gene expression profiling and meta-analysis. Nucleic Acids Res (2019) 47(W1):W234–W41. doi: 10.1093/nar/gkz240

PubMed Abstract | CrossRef Full Text | Google Scholar

28. O'Neill TW, McCabe PS, McBeth J. Update on the epidemiology, risk factors and disease outcomes of osteoarthritis. Best Pract Res Clin Rheumatol (2018) 32(2):312–26. doi: 10.1016/j.berh.2018.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Bortoluzzi A, Furini F, Scire CA. Osteoarthritis and its management - epidemiology, nutritional aspects and environmental factors. Autoimmun Rev (2018) 17(11):1097–104. doi: 10.1016/j.autrev.2018.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Hwang HS, Kim HA. Chondrocyte apoptosis in the pathogenesis of osteoarthritis. Int J Mol Sci (2015) 16(11):26035–54. doi: 10.3390/ijms161125943

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Weng PW, Pikatan NW, Setiawan SA, Yadav VK, Fong IH, Hsu CH, et al. Role of gdf15/mapk14 axis in chondrocyte senescence as a novel senomorphic agent in osteoarthritis. Int J Mol Sci (2022) 23(13):7043. doi: 10.3390/ijms23137043

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Liu D, Liu W, Jiang L, Dong S, Ma W, Wang S, et al. Silencing of tlr7 protects against lipopolysaccharide-induced chondrocyte apoptosis and injury by blocking the P21-mediated jak2/stat3 pathway. Am J Transl Res (2021) 13(12):13555–66.

PubMed Abstract | Google Scholar

34. Hosseinzadeh A, Kamrava SK, Joghataei MT, Darabi R, Shakeri-Zadeh A, Shahriari M, et al. Apoptosis signaling pathways in osteoarthritis and possible protective role of melatonin. J Pineal Res (2016) 61(4):411–25. doi: 10.1111/jpi.12362

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Snead MP, Yates JR. Clinical and molecular genetics of stickler syndrome. J Med Genet (1999) 36(5):353–9. doi: 10.1136/jmg.36.5.353

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Kannu P, Bateman JF, Randle S, Cowie S, du Sart D, McGrath S, et al. Premature arthritis is a distinct type ii collagen phenotype. Arthritis Rheum (2010) 62(5):1421–30. doi: 10.1002/art.27354

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Ruan G, Ying Y, Lu S, Zhu Z, Chen S, Zeng M, et al. The effect of systemic iron status on osteoarthritis: A mendelian randomization study. Front Genet (2023) 14:1122955. doi: 10.3389/fgene.2023.1122955

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Valdes AM, Spector TD. Genetic epidemiology of hip and knee osteoarthritis. Nat Rev Rheumatol (2011) 7(1):23–32. doi: 10.1038/nrrheum.2010.191

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Boer CG, Hatzikotoulas K, Southam L, Stefansdottir L, Zhang Y, Coutinho de Almeida R, et al. Deciphering osteoarthritis genetics across 826,690 individuals from 9 populations. Cell (2021) 184(18):4784–818.e17. doi: 10.1016/j.cell.2021.07.038

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Bi J, Li Q, Yang Z, Cai L, Lv T, Yang X, et al. Cxcl2 impairs functions of bone marrow mesenchymal stem cells and can serve as a serum marker in high-fat diet-fed rats. Front Cell Dev Biol (2021) 9:687942. doi: 10.3389/fcell.2021.687942

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Manoharan P, Song T, Radzyukevich TL, Sadayappan S, Lingrel JB, Heiny JA. Klf2 in myeloid lineage cells regulates the innate immune response during skeletal muscle injury and regeneration. iScience (2019) 17:334–46. doi: 10.1016/j.isci.2019.07.009

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Wang Z, He L, Li W, Xu C, Zhang J, Wang D, et al. Gdf15 induces immunosuppression via cd48 on regulatory T cells in hepatocellular carcinoma. J Immunother Cancer (2021) 9(9):e002787. doi: 10.1136/jitc-2021-002787

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Tung CB, Li CY, Lin HY. Multi-omics reveal the immunological role and the theragnostic value of mir-216a/gdf15 axis in human colon adenocarcinoma. Int J Mol Sci (2021) 22(24):13636. doi: 10.3390/ijms222413636

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Yang SL, Tan HX, Lai ZZ, Peng HY, Yang HL, Fu Q, et al. An active glutamine/alpha-ketoglutarate/hif-1alpha axis prevents pregnancy loss by triggering decidual igf1(+)Gdf15(+)Nk cell differentiation. Cell Mol Life Sci (2022) 79(12):611. doi: 10.1007/s00018-022-04639-x

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Husaini Y, Tsai VW, Manandhar R, Zhang HP, Lee-Ng KKM, Lebhar H, et al. Growth differentiation factor-15 slows the growth of murine prostate cancer by stimulating tumor immunity. PloS One (2020) 15(6):e0233846. doi: 10.1371/journal.pone.0233846

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Pattappa G, Schewior R, Hofmeister I, Seja J, Zellner J, Johnstone B, et al. Physioxia has a beneficial effect on cartilage matrix production in interleukin-1 beta-inhibited mesenchymal stem cell chondrogenesis. Cells (2019) 8(8):936. doi: 10.3390/cells8080936

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Jenei-Lanzl Z, Meurer A, Zaucke F. Interleukin-1beta signaling in osteoarthritis - chondrocytes in focus. Cell Signal (2019) 53:212–23. doi: 10.1016/j.cellsig.2018.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Ciucci T, Ibanez L, Boucoiran A, Birgy-Barelli E, Pene J, Abou-Ezzi G, et al. Bone marrow th17 tnfalpha cells induce osteoclast differentiation, and link bone destruction to ibd. Gut (2015) 64(7):1072–81. doi: 10.1136/gutjnl-2014-306947

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Ishii H, Tanaka H, Katoh K, Nakamura H, Nagashima M, Yoshino S. Characterization of infiltrating T cells and th1/th2-type cytokines in the synovium of patients with osteoarthritis. Osteoarthritis Cartil (2002) 10(4):277–81. doi: 10.1053/joca.2001.0509

CrossRef Full Text | Google Scholar

50. Won HY, Lee JA, Park ZS, Song JS, Kim HY, Jang SM, et al. Prominent bone loss mediated by rankl and il-17 produced by cd4+ T cells in tallyho/jngj mice. PloS One (2011) 6(3):e18168. doi: 10.1371/journal.pone.0018168

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Saravia J, Chapman NM, Chi H. Helper T cell differentiation. Cell Mol Immunol (2019) 16(7):634–43. doi: 10.1038/s41423-019-0220-6

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Charo IF, Ransohoff RM. The many roles of chemokines and chemokine receptors in inflammation. N Engl J Med (2006) 354(6):610–21. doi: 10.1056/NEJMra052723

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Walsh MC, Choi Y. Regulation of T cell-associated tissues and T cell activation by rankl-rank-opg. J Bone Miner Metab (2021) 39(1):54–63. doi: 10.1007/s00774-020-01178-y

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Croes M, Oner FC, van Neerven D, Sabir E, Kruyt MC, Blokhuis TJ, et al. Proinflammatory T cells and il-17 stimulate osteoblast differentiation. Bone (2016) 84:262–70. doi: 10.1016/j.bone.2016.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Hu Y, Ek-Rylander B, Wendel M, Andersson G. Reciprocal effects of interferon-gamma and il-4 on differentiation to osteoclast-like cells by rankl or lps. Oral Dis (2014) 20(7):682–92. doi: 10.1111/odi.12189

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Gerritsen B, Pandit A. The memory of a killer T cell: models of cd8(+) T cell differentiation. Immunol Cell Biol (2016) 94(3):236–41. doi: 10.1038/icb.2015.118

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Toni R, Di Conza G, Barbaro F, Zini N, Consolini E, Dallatana D, et al. Microtopography of immune cells in osteoporosis and bone lesions by endocrine disruptors. Front Immunol (2020) 11:1737. doi: 10.3389/fimmu.2020.01737

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Metcalfe DD. Mast cells and mastocytosis. Blood (2008) 112(4):946–56. doi: 10.1182/blood-2007-11-078097

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Al-Hakami A, Alqhatani SQ, Shaik S, Jalfan SM, Dhammam MSA, Asiri W, et al. Cytokine physiognomies of mscs from varied sources confirm the regenerative commitment post-coculture with activated neutrophils. J Cell Physiol (2020) 235(11):8691–701. doi: 10.1002/jcp.29713

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Harrell CR, Markovic BS, Fellabaum C, Arsenijevic A, Volarevic V. Mesenchymal stem cell-based therapy of osteoarthritis: current knowledge and future perspectives. BioMed Pharmacother (2019) 109:2318–26. doi: 10.1016/j.biopha.2018.11.099

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Chen K, Jiao Y, Liu L, Huang M, He C, He W, et al. Communications between bone marrow macrophages and bone cells in bone remodeling. Front Cell Dev Biol (2020) 8:598263. doi: 10.3389/fcell.2020.598263

PubMed Abstract | CrossRef Full Text | Google Scholar