Bioinformatics-Led Discovery of Osteoarthritis Biomarkers and Inflammatory Infiltrates

The molecular mechanisms of osteoarthritis, the most common chronic disease, remain unexplained. This study aimed to use bioinformatic methods to identify the key biomarkers and immune infiltration in osteoarthritis. Gene expression profiles (GSE55235, GSE55457, GSE77298, and GSE82107) were selected from the Gene Expression Omnibus database. A protein-protein interaction network was created, and functional enrichment analysis and genomic enrichment analysis were performed using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genome (KEGG) databases. Immune cell infiltration between osteoarthritic tissues and control tissues was analyzed using the CIBERSORT method. Identify immune patterns using the ConsensusClusterPlus package in R software using a consistent clustering approach. Molecular biological investigations were performed to discover the important genes in cartilage cells. A total of 105 differentially expressed genes were identified. Differentially expressed genes were enriched in immunological response, chemokine-mediated signaling pathway, and inflammatory response revealed by the analysis of GO and KEGG databases. Two distinct immune patterns (ClusterA and ClusterB) were identified using the ConsensusClusterPlus. Cluster A patients had significantly lower resting dendritic cells, M2 macrophages, resting mast cells, activated natural killer cells and regulatory T cells than Cluster B patients. The expression levels of TCA1, TLR7, MMP9, CXCL10, CXCL13, HLA-DRA, and ADIPOQSPP1 were significantly higher in the IL-1β-induced group than in the osteoarthritis group in an in vitro qPCR experiment. Explaining the differences in immune infiltration between osteoarthritic tissues and normal tissues will contribute to the understanding of the development of osteoarthritis.

Introduction

Osteoarthritis is a degenerative joint disease that affects the elderly the most (1). It is the major cause of disability and is characterized by cartilage destruction and bone fragmentation in joints (2, 3). Prior joint injury, abnormal joint or limb development, genetic background, and having a job that requires heavy lifting are among the risk factors for osteoarthritis caused by mechanical stress on the joints and low levels of inflammation (4). osteoarthritis is the most common chronic joint disease, which cannot be prevented, and its prevalence rises with age (5–7). The treatment costs for osteoarthritis have created a financial burden for patients and society, causing a serious threat to human health. The disease’s etiology and pathogenesis, however, are still unknown. To develop effective treatment strategies, it is necessary to investigate the underlying mechanisms of osteoarthritis development.

Animal models (7), tissue models (8), gene expression (9), and systems biology approach (10) have been used to investigate the molecular mechanisms of osteoarthritis. Many researchers have focused on the identification of disease-associated proteins for osteoarthritis, such as type 2 collagen(COL2A1) and aggregated proteoglycans(ACAN) (11), nucleotide-binding and oligomerized structural domain receptor inflammatory vesicles containing protein 3 (NLRP3) (12), and matrix metalloproteinase (MMP)-13 (13). The molecular basis of Osteoarthritis pathology, however, remains unexplained.

The involvement of differentially expressed genes (DEGs) and biofunctional pathways, which are analyzed using microarray and bioinformatic technologies, in the development of osteoarthritis microarray and bioinformatics has extended the scope of the previous genome-based studies to screening for genetic alterations (10) Bioinformatics is an interdisciplinary method used to study the molecular mechanisms of diseases (14). Improved knowledge of the molecular networks and genes involved in those networks will enhance the overall understanding of the pathology of osteoarthritis and the molecular networks and genes involved. High false positive rates were observed in the study that used only one microarray platform and a small number of samples, which may have contributed to the inconsistent results (15). Therefore, further research is required to identify new therapeutic targets and diagnostic biomarkers that are more reliable to overcome the inconsistencies observed in previous studies.

In this study, gene expression profiles and microarray experiments were used to analyze the genes that were significantly different between patients with osteoarthritis and control samples. The identification of key biomarkers of unusually expressed genes and immune infiltration will improve our understanding of the molecular mechanisms of osteoarthritis at the systems biology level.

Materials and Methods

Microarray Data Source

The datasets GSE55235 (16), GSE55457 (16), GSE77298 (17), and GSE82107 (18) were downloaded from the Gene Expression Omnibus (GEO) database. Microarray data of GSE55235 contain 20 arthritis samples and 10 control samples, those of GSE55457 contain 23 arthritis samples and 10 control samples, those of GSE77298 contain 16 arthritis samples and 10 control samples, and those of GSE82107 contain 16 arthritis samples and 10 control samples. The GSE55235 and GSE55457 datasets were sequenced on GPL96, GSE77298 and GSE82107 on GPL570, and GSE77298 and GSE82107 on GPL570, all with the human body as the origin. (Table 1) To create the integrated GEO dataset, the four datasets mentioned above were de-batched using the R package sva (19) to contain 69 arthritis samples and 34 control samples.

TABLE 1

Table 1 Descriptive statistics.

The ImmPort (20), GenegCards (21), and MSigDB (22) databases were used to obtain 1509, 16,664, and 21,341 immune-related genes, respectively. Furthermore, 1264 immune-related genes were derived from the intersection of the aforementioned three gene sets.

Identification of Arthritis-Related Immune-Differentially Expressed Genes

Differential gene analysis was performed using the R package “limma” (23) to determine the differential genes between diseased and control samples in the integrated dataset to investigate the impact of gene expression levels of immune-related genes on arthritis. The thresholds for differential genes were set at the absolute value of log2 fold change |log2FC| >1 and adj. P <0.05, indicating DEGs with upregulated expression. Similarly, log2FC <-1and adj. P <0.05 indicates DEGs with downregulated expression. Volcano plots and heat maps were used to display the results of differential gene expression. The differentially expressed immune genes were identified by intersecting DEGs and immune genes.

Forest Model and Nomogram Model Construction

The model was trained using the least absolute shrinkage and selection operator (LASSO) technique to predict the likelihood of arthritis. The candidate differentially expressed immune genes were included in the model and analyzed using the LASSO algorithm to obtain the characteristic genes associated with arthritis. The risk score formula (Eq. 1) was established using the forest model to predict the likelihood of arthritis.

$\begin{array}{ll}\text{Risk\hspace{0.17em}Score}=\left(\text{exp}-\text{Gene}1*\text{coef}-\text{Gene}1\right)+\left(\text{exp}-\text{Gene}2*\text{coef}-\text{Gene}2\right)+\dots +\left(\text{exp}-\text{Gene}*\text{coef}-\text{Gene}\right)& (1)\end{array}$

The prevalence of arthritis in patients was then predicted using columnar line graph models based on the chosen candidate Imm-DEGs. The expected values were compared to the standard value using calibration curves. To determine whether the model-based decisions were beneficial to patients, decision curve analysis was performed and clinical impact curves were drawn.

Identification of Molecular Subtypes Based on Important Immunomodulators

Consistency clustering is a resampling-based approach for identifying each member and their subgroup number, as well as validating the cluster. To discover various immunological patterns based on significant Imm-DEGs, the ConsensusClusterPlus package (24) in R was employed.

Determination of DEGs in Various Immunological Patterns

To investigate the effect of diverse immunological patterns on arthritis, differentially expressed genes were screened for significant differential genes using the R package limma on immune pattern subgroups in the integrated dataset. Differentially expressed genes with |log2FC| >1 and adj. P <0.05 have upregulated expression, while those with log₂FC=-1 and adj. P <0.05 has downregulated expression. The results of the differential gene expression were obtained using volcano and heat maps.

Assessment of Biological Variables Among Various Immunological Patterns

Gene Ontology (GO) analysis is a standard method to conduct large-scale functional enrichment studies, covering biological processes, molecular functions, and cellular components (24). Kyoto Encyclopedia of Genes and Genome (KEGG) is a widely used database for storing information on genomes, biological pathways, diseases, and pharmaceuticals (25). GO annotation analysis and KEGG pathway enrichment analysis of differentially expressed genes were performed using R’s clusterProfiler package (26), and a threshold value of <0.05 for false discovery rate was considered statistically significant.

To investigate the differences in biological processes between different subgroups, based on a dataset of gene expression profiles from arthritis patients, we performed gene set enrichment analysis (GSEA). GSEA is a computational method used to analyze whether a particular gene set is statistically significant and consistent with the differences between two biological states (22). GSEA is commonly used to estimate the changes in the pathways and biological activities in the samples of expression datasets.The gene sets “c2.cp.kegg.v7.4.symbols.gmt” and “c5.go.v7.2.symbols.gmt” were retrieved from the MSigDB database for use in the GSEA. P-value <0.05 was considered statistically significant.

Protein-Protein Interaction Network

Protein-protein interaction (PPI) networks are created for individual proteins interacting with one other to participate in numerous activities of life processes, such as biological signaling, control of gene expression, energy and material metabolism, and cell cycle regulation. Systematic analysis of the interactions of a large number of proteins in biological systems is important to understand how proteins work in biological systems. Systematic analysis is also useful to understand the response mechanisms of biological signaling and energy-matter metabolism in specific physiological states (disease) and functional connections between proteins. The STRING database (27) is widely used for searching for known proteins and predicting relationships between proteins. Using the STRING database, we constructed a PPI network linked with differentially expressed genes, Imm-DEGs, and prospective differentially expressed immune genes. The PPI network model was displayed using Cytoscape (v3.7.2) (28), and the genes in the network were functionally annotated using closeGO (29).

Identification and Correlation of Disease Immune Infiltrate Cells

The immune microenvironment generally consists of immune cells, inflammatory cells, fibroblasts, mesenchymal samples, different cytokines, and chemokines. Immune cell infiltration analysis has an important guiding role in the prediction of the disease course and response to therapy. The single sample gene set enrichment analysis (ssGSEA) algorithm, an extension of the GSEA approach, was employed to quantify the abundance of 28 immune cell types in individuals with different immunological patterns (30). The CIBERSORT algorithm, which can perform linear support vector regression to deconvolute gene expression profiles, was used to estimate the number of immune cells in samples using RNA-sequencing data (31). We calculated 22 immune cell types in patients with distinct immunological patterns in the dataset using the CIBERSORT algorithm (31) in R software and showed the immune cell composition of patients with varied immune patterns using box plots. Differences in immune cell proportions were evaluated using the Wilcoxon rank-sum test. P <0.05 was considered statistically significant.

qRT-PCR Validation of the Hub Genes

Interleukin-1β (IL-1β) can induce cartilage degradation by promoting the expression of matrix metalloproteinases (MMPs) in chondrocytes (32) and is widely used in the inflammatory induction model of chondrocytes in osteoarthritis (33) . We used IL-1β to stimulate chondrocytes of normal and osteoarthritis person to simulate this microenvironment of inflammation and used qRT-PCR to verify the expression level of hub genes. The chondrocytes from normal people (CP-H107, Procell, wuhan, Hubei, China) and patients with osteoarthritis (402OAK-05a, Haoge Biotechnology Co., Ltd, Shanghai, China) and were cultured in DMEM/F12 medium (SH30126.01, HyClone Technologies, Logan, USA) containing 10% fetal bovine serum (FBS, 10099, Thermo Fisher Scientific, Massachusetts, USA) in 5% CO₂ at 37°C and divided into NC (negative control) group and IL-1β group. Cells in the NC group were added to the normal medium while cells in the IL-1β group were treated with IL-1β (40 ng/mL) culture medium. The mRNA relative expression of TCA1, TLR7, MMP9, CXCL10, CXCL13, HLA-DRA, ADIPOQ, LEP, and SPP1 were detected after 48 h of incubation. The total RNA of osteoarthritis was extracted and synthesized into cDNA according to the manufacturer’s protocol (Promega Biotech Co., Ltd, Beijing, China). qRT-PCR was performed on a LightCycler 96 (Roche Life Sciences, Switzerland, Basel) using Real-Time PCR Mix (Vazyme Biotech, Nanjing, Jiangsu Province, China). Gene expression relative to GAPDH expression was assessed using the 2^-ΔΔCt method. Independent experiments were conducted in triplicate. (The sequence fragments of RNAs are shown in Table 2).

TABLE 2

Table 2 PCR primers.

The experiment was divided into 2 groups: Interleukin-1β (IL-1β) group cells were added to a drug-containing medium of 40 ng/mL; osteoarthritis group cells were added to the normal medium. The mRNA relative expression of TCA1, TLR7, MMP9, CXCL10, CXCL13, HLA-DRA, ADIPOQ, LEP, and SPP1 were detected after 48 h of incubation. The primers were designed using the DNAMAN software and synthesized by Sangon Biotech (Shanghai, China), with primer sequences shown in Table 2. The cellular RNA was extracted using TRIzol (Invitrogen #15596-026). cDNA was synthesized using PrimeScript™ RT reagent Kit with gDNA Eraser (Takara#RR047A) and SYBR Green qPCR Mix (Beyotime#D7260). The 7500 Real-Time Polymerase Chain Reaction (RT-PCR) System was used to complete amplification in 40 cycles. The PCR data were processed using GAPDH as an internal reference, and the relative expression in the samples was calculated using the AGCT method.

Statistical Analysis

All data processing and analysis were performed in RStudio (version 4.1.1). To compare two groups of continuous variables, statistical significance of normally distributed variables was calculated using an independent Student’s t-test, and differences between non-normally distributed variables were calculated using the Mann-Whitney U-test (i.e., Wilcoxon rank-sum test). The chi-square test or Fisher’s exact test was carried out to analyze the statistical significance between two sets of categorical variables. Correlation coefficients between different genes were estimated via Pearson correlation analysis. All statistical P values were two-sided, and p <0.05 was considered statistically significant.

Results

Expression of Immune-Related Genes in Arthritis Patients

The bioinformatics analysis of this study is carried out according to Figure 1. The batch effects were removed from the GEO dataset to obtain the integrated dataset (Figure 2), which includes 69 arthritis samples and 34 control samples. Differential study of arthritic samples and control samples revealed 105 DEGs; 83 of them were upregulated and 22 were downregulated (Figures 3A, B). Figure 3C shows the intersection of immune-related genes from three datasets and the intersection of DEGs on immune-related genes. These intersections yielded 28 Imm-DEGs (Figure 3D), 26 of which had an upregulated expression pattern (Figure 3E) and 2 had a downregulated one (Figure 3F).

FIGURE 1

Figure 1 Flow chart.

FIGURE 2

Figure 2 GEO data de-batching. (A) Gene expression level statistics of the dataset before de-batching. (B) Gene expression level statistics of the integrated dataset after de-batching. (C) Principal component analysis (PCA) between datasets before de-batching. (D) PCA between integrated datasets after de-batching.

FIGURE 3

Figure 3 Immune-related genes and differentially expressed immune genes (Imm-DEGs). (A) Arthritis-related differentially expressed genes (DEGs) volcano plot with log2FoldChange in the horizontal coordinate and -log₁₀(Adjust P-value) in the vertical coordinate. Red nodes indicate upregulated DEGs, blue nodes indicate downregulated DEGs, and gray nodes indicate genes that are not significantly differentially expressed. (B) Heat map of arthritis-related DEG expression levels: pink indicates disease samples, green indicates normal control samples, red indicates high gene expression, and blue indicates low gene expression. (C) Immune gene Venn diagram: three colors represent three different data sources. (D) Immune gene versus DEG Venn diagram: pink represents immune genes, green represents DEGs (E) Venn diagram of immune genes and upregulated DEGs: pink represents immune genes, green represents upregulated DEGs. (F) Venn diagram of immune genes and downregulated DEGs: pink represents immune genes, green represents downregulated DEGs.

To analyze the overall expression of Imm-DEGs, heat maps (Figure 4A) and histograms of the expression levels of Imm-DEGs in arthritis samples and control samples were drawn. Most genes were expressed at higher levels in arthritis samples than in normal samples (Figure 4C), and the positions of Imm-DEGs were annotated on human chromosomes (Figure 4B).

FIGURE 4

Figure 4 Expression levels of Imm-DEGs in arthritis. (A) Heat map of overall expression of Imm-DEGs in arthritis patients: green for control samples, pink for disease samples, red for high expression, blue for low expression. (B) Chromosome distribution of immune-related genes in arthritis patients. (C) Overall expression histogram of immune-related genes in arthritis patients: green for control samples, pink for disease samples, horizontal axis indicates genes, vertical axis indicates gene expression levels. (***P < 0.001).

Construction of Risk Model

The LASSO algorithm was used to identify 16 characteristic genes out of 28 Imm-DEGs with a great effect on arthritis (Figures 5A, B). Based on the coefficients of 16 characteristic genes (Figure 5C), the gene expressions were multiplied by the corresponding coefficients and summed to derive the arthritogenic score (Figure 5D). Similarly, the receiver operating characteristic (ROC) curves of 16 gene signatures were investigated to predict arthritis, and the findings revealed the predictive efficacy of all 16 gene signatures (Figure 5E).

FIGURE 5

Figure 5 Construction of arthritis model. (A, B) Screening of gene signatures from Imm-DEGs using the LASSO algorithm. (C) Forest plot of gene signatures in arthritis patients. (D) Receiver operating characteristic (ROC) curve of predicted risk scores in arthritis diagnosis. (E) ROC curve of 16 gene signatures in arthritis diagnosis.

A line plot model was created based on the patient’s projected risk score and 16 trait genes to predict the prevalence of arthritis in patients (Figures 6A, B). The calibration curves revealed that the line graph model predictions were nearly identical to those of the ideal model (Figure 6C) and that the single predicted risk score in the decision curve analysis or the composite genetic model is better than that in the random model. These results imply that decision-making based on the line graph model could be beneficial for arthritis patients (Figure 6D).

FIGURE 6

Figure 6 Line graphs (nomogram plots). (A) Nomogram of predicted risk scores in the diagnosis of arthritis patients. (B) Nomogram of 16 trait genes in the diagnosis of arthritis patients. (C) Nomo model evaluation, where the diagnostic model is in better agreement with the ideal model. (D) Model evaluation curves: gray indicates the follow-on diagnosis, green indicates the complex diagnostic model with 17 trait genes, and pink indicates the simple diagnostic model with predicted risk score. (*P<0.05, ***P<0.001).

Correlations between gene expression levels and functional correlations between 16 gene signatures were examined. Coefficients of functional correlations between IGKV1-17, IGLC1, and other genes reached 0.9 (Figure 7B). In all samples, the correlation coefficient between IGKV1-17 and IGLC1 was 0.95; ADIPOQ had a positive correlation with LEP but a negative correlation with all other genes except ADIPOQ, and LEP had a negative correlation with all other genes except ADIPOQ (Figure 7A). In arthritic samples, the correlation coefficient between IGKV1-17 and IGLC1 was 0.94; ADIPOQ showed a positive correlation with LEP and SPP1 and a negative correlation with all other genes except ADIPOQ. Similarly, LEP showed a negative correlation with all other genes except ADIPOQ (Figure 7C). OGN gene also showed a negative correlation with most other characteristic genes. In control samples, the correlation coefficient between IGKV1-17 and IGLC1 was 0.89 (Figure 7D).

FIGURE 7

Figure 7 (A) Correlation analysis of 16 genes in all samples: * represents the significance of the correlation, and the number represents the correlation level. (B): Functional correlation analysis of 16 trait genes, the horizontal axis indicates the correlation size and the vertical axis indicates the trait genes. (C, D) Correlation analysis of 16 trait genes in disease and normal samples: * represents the significance of correlation, and the number represents the correlation level. (*P<0.05, **P<0.01, ***P<0.001).

Distinct Immunological Patterns of Gene Signatures

Two immunological patterns (clusterA and clusterB) were established utilizing the ConsensusClusterPlus package in R software and a consistent clustering approach based on 16 signature genes. There were 34 samples in Cluster A and 35 samples in Cluster B. Heat maps for all differentially expressed immune genes were then created to show the considerable differences in immune gene expression between the two groupings (Figure 8A). Expression levels of ADIPOQ, LEP, OGN, and TAC1 were significantly lower in cluster A than that in cluster B, whereas HLA-DMA, HLA-DMB, HLA-DPA1, HLA-DQA1, HLA-DQB1 HLA-DRA, TLR7, CCL18, CXCL10, MMP9, BLNK, IGHM, IGKV1-17, IGLC1, IGLV1-44, CXCL13, CXCL6, CXCL9, SPP1, IL10RA, SDC1, TNFRSF17, and TRBC1 expression levels were significantly higher in cluster A than that in cluster B (Figure 8B). Similarly, the ROC curves of the 16 gene signatures individually predicted for both categories were evaluated, and the findings showed that all 16 gene signatures had good classification efficacy (Figure 8C).

FIGURE 8

Figure 8 Consistency clustering of gene signatures in arthritis patients. (A) Heat map of expression levels of 28 Imm-DEGs in two clusters: pink indicates cluster A, green indicates cluster B, red indicates high expression, and blue indicates low expression. (B) Expression levels of 28 Imm-DEGs in two clusters: pink indicates cluster A, green indicates cluster B, the horizontal axis is the Imm-DEGs, and the vertical axis is the gene expression level. (C) ROC curves of 16 characterized genes independently distinguished between cluster A and cluster B. (nsP≥0.05, *P<0.05, **P<0.01, ***P<0.001).

PPI Network of Immune Genes

To explore the relationship between differentially expressed immune genes, we extracted the PPI network of DEGs, Imm-DEGs, and gene signatures. As visualized in Cytoscape, the PPI network of DEGs had 211 pairing interactions and 75 genes; MMP9 was strongly correlated with 19 DEGs, whereas CXCL10 was linked to 16 DEGs (Figure 9A). Similarly, the PPI network of Imm-DEGs comprised 58 reciprocal pairs and 24 genes; CXCL10 and MMP9 were closely linked to 10 differentially expressed immune genes, whereas CXCL13 and HLA-DRA were both linked to 8 differentially expressed immune genes (Figure 9C). The PPI network of the gene signatures contained 25 interaction pairs and 15 genes, where CXCL10 and MMP9 were closely linked to 10 Imm-DEGs, whereas CXCL13 and HLA-DRA were both linked to 8 Imm-DEGs (Figure 9E). To verify the functions of genes in three PPI networks, functional enrichment analysis was performed using ClueGO. The results revealed that genes in the PPI network were enriched in the pathways of adipocytokine signaling, glycolysis/gluconeogenesis, toll-like receptor signaling, viral protein interaction with cytokine and cytokine receptor IL-17 signaling, and rheumatoid arthritis (Figure 9B). Gene enrichment in the PPI network for differentially expressed immune genes was involved in the intestinal immune network for IgA production, viral protein interaction with cytokine and cytokine receptors (Figure 9D).

FIGURE 9

Figure 9 Immune-related gene protein-protein interaction (PPI) network. (A) Differentially expressed gene PPI network: yellow nodes indicate Imm-DEGs. (B) Results of gene enrichment analysis in DEG PPI network. (C) Differentially expressed immune gene PPI network: yellow nodes indicate gene signatures. (D) Results of gene enrichment analysis in Imm-DEG PPI network. (E) Results of gene signatures.

Differential Analysis of Two Different Immune Patterns

To analyze the differences between the two immune patterns, 168 DEGs were obtained between the patterns: cluster A and cluster B, including 93 DEGs upregulated and 75 DEGs downregulated in cluster A (Figure 10A). The heat maps showed that these DEGs could distinguish between the two different immune patterns (Figure 10B).

FIGURE 10

Figure 10 Differential analysis between two different immune patterns. (A) Horizontal coordinate is log2FoldChange; vertical coordinate is -log₁₀ (Adjust P-value); red nodes indicate upregulated DEGs; blue nodes indicate downregulated DEGs; gray nodes indicate genes that are not significantly differentially expressed. (B) Heat map of expression levels of DEGs in two clusters: pink indicates cluster A; blue indicates cluster B; red indicates high expression; blue indicates low expression.

Subsequently, we analyzed the role of DEGs between the two immune modalities in the biologically relevant functions of patients. First, the DEGs were functionally annotated (Figure 11A), Second, these DEGs were enriched in cytokine-cytokine receptor interaction, chemokine signaling, malaria, and tumor necrosis factor (TNF) signaling, and NOD-like receptor signaling pathways upon KEGG pathway analysis (Figure 11B).

FIGURE 11

Figure 11 Functional analysis between two different immune patterns. (A) Gene Ontology (GO) functional enrichment analysis: vertical coordinate is the significance of the enrichment result; horizontal coordinate is the Z-score; node color indicates BP, CC, MF; node size indicates the number of genes contained in the current GO term. (B) Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis results: node color indicates gene expression level; quadrilateral color indicates Z-score. BP, Biological Process; MF, molecular function; CC, cellular component.

Next, we performed GSEA on all genes between two immune modalities, which were significantly different in biological processes such as activation of immune response in cluster A, adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains, alpha-beta T cell activation, and antigen binding. On the contrary, biological processes such as primary alcohol metabolic processes, regulation of insulin-like growth factor receptor signaling pathway, lipid oxidation, and retinol metabolism were inhibited (Figure 12A). Chemokine signaling pathway, cytokine receptor interaction, graft versus host disease, leishmanial infection, and NK cell-mediated cytotoxicity were activated, whereas butanoate metabolism, metabolism of xenobiotics by cytochrome p450, fatty acid metabolism, valine leucine and isoleucine degradation, and drug metabolism cytochrome p450 pathways were inhibited (Figure 12B).

FIGURE 12

Figure 12 Gene Set Enrichment Analysis (GSEA) between two different immune models. (A) GSEA-GO analysis: horizontal coordinate is the gene ratio; vertical coordinate shows the GO terms; color indicates -log₁₀ (p-value); node size indicates the number of genes enriched in GO terms. (B) GSEA-KEGG analysis: horizontal coordinate is the gene ratio; vertical coordinate shows the GO terms; node size indicates the number of genes enriched in the pathway; node color indicates -log₁₀ (p-value).

Differences in Immune Characteristics Between Two Models

The CIBERSORT algorithm was used to assess the level of immune cell infiltration between two different immune modalities. CIBERSORT analysis revealed that patients in cluster A had significantly lower levels of resting dendritic cells, M2 macrophages, resting mast cells, activated NK cells, and regulatory T cells than those in cluster B (Figure 13A). However, the levels of M1 macrophages, activated mast cells, plasma cells, T follicular helper cells, and gamma delta T cells were significantly higher than those in cluster B (Figure 13B).

FIGURE 13

Figure 13 Immune characteristics between two different immune patterns. (A) Immune cell content stacking plot between cluster A and cluster B: different colors indicate different immune cells; the horizontal axis is the patient ID. (B) Immune cell content histogram: the horizontal axis indicates 22 immune cells; the vertical axis indicates cell content; pink indicates cluster A samples; green indicates cluster B samples. (C) Immune cell content histogram: the horizontal axis indicates 22 immune cells; the vertical axis indicates the cell content; pink indicates cluster A samples; green indicates cluster B samples. Correlations between 12 immune cells that significantly differed between cluster A and cluster B patients in the limma algorithm: correlations in cluster A patient species (C); correlations in cluster B patient species (D); red indicates negative correlations; blue indicates positive correlations. Correlations in all arthritis patients. Correlation analysis of 16 gene signatures and immune cell content in all arthritis patients (E), cluster A arthritis patients (F) and cluster B arthritis patients (G), horizontal axis indicates immune cells, vertical axis indicates 16 gene signatures, node color indicates the correlation size, and node size indicates the significance level. (?:no data, nsP≥0.05, **P<0.01, ***P<0.001).

The correlation between the immune cell contents of patients in clusters A and B was also calculated. The patients in cluster A showed a higher proportion of T follicular helper cells than that naïve B cells. Gamma delta T cells were in higher proportion than M1 macrophages. Furthermore, T follicular helper cells and activated memory CD4 T cells showed a significant positive correlation, whereas activated and resting mast cells, M2 macrophages, and naive B cells showed a significant negative correlation (p < 0.05, Figure 13C). Similarly, in cluster B patients, gamma delta T cells, M1 macrophages, T follicular helper cells, gamma delta T cells, and activated mast cells showed a positive correlation with regulatory T cells, whereas regulatory T cells showed a negative correlation with M1 macrophages, regulatory T cells, and gamma delta T cells (p <0.05, Figure 13D).

The correlations between 16 gene signatures and immune cell types were also calculated for all arthritis patients (Figure 13E), cluster A (Figure 13F), and cluster B arthritis patients (Figure 13G). TAC1 was not correlated with the immune cell content in cluster B arthritis patients, and resting memory CD4 T cells and activated NK cells were not correlated with the gene signatures (Figure 13G).

qPCR Validation of Data

To verify the bioinformatics results, qPCR experiments were conducted. The results revealed that the mRNA expression levels of TCA1, TLR7, MMP9, CXCL10, CXCL13, HLA-DRA, and ADIPOQSPP1 were significantly higher in the IL-1β-induced group, with the most significant difference in the CXCL10 expression. The difference in the LEP expression between the two groups was not statistically significant. This indicates that the results of data mining are reliable and have potential research value (Figure 14).

FIGURE 14

Figure 14 qPCR validation. After IL-1β-induced, the mRNA expression levels of TCA1, TLR7, MMP9, CXCL10, CXCL13, HLA-DRA, ADIPOQ and SPP1 were significantly higher both in the normal person (A) and osteoarthritis (B). (nsP≥0.05, *P<0.05).

Discussion

Microarray technology and high-throughput technology, the main approaches to exploring the expression levels of genes, have improved the understanding of the intrinsic molecular mechanisms associated with complex disorders. Osteoarthritis is a degenerative joint condition that commonly affects the hands and weight-bearing joints (34). It is the most common joint disease in the world (35), affecting approximately 10–15% of the youth and 60% of the elderly (36). The molecular mechanism of its immune milieu is poorly understood; therefore, finding a novel target for early diagnosis, treatment, and prognosis of the immune microenvironment of osteoarthritis has important clinical benefits. All arthritis and control samples in the present study comprised 105 differently expressed genes; 28 of them were immune-related genes that overlapped with differentially expressed genes. Inflammatory response, immunological response, and the chemokine-mediated signaling pathway were the most highly enriched categories in the GO analysis of DEGs. Differentially expressed genes were found to be involved in cytokine-cytokine receptor interactions, chemokine signaling pathways, and malaria, as revealed by the KEGG analysis. Cluster A patients had significantly lower levels of resting dendritic cells, M2 macrophages, resting mast cells, and activated NK cells than cluster B patients did, whereas M1 macrophages, activated mast cells, plasma cells, and T follicular helper cells in cluster A were significantly higher than those in cluster B patients.

In the PPI reciprocal network analysis, CXCL13, CXCL9, SPP1, and CXCL10 were linked to Imm-DEGs. Osteoarthritic osteoblasts directly regulate proliferation and type I collagen expression by CXCL13 chemokines (37). As reported in rheumatoid arthritis, proinflammatory cytokines enhance secondary CXCL13 production by reactivating CXCL13-producing CD4 T cells (38). Similarly, serum CXCL13 correlates with synovial CXCL13 measured at the joints, suggesting that synovitis is an important source of circulating CXCL13 (39). Past evidence showed the linkage between tissue bone bridging protein (SPP1) and osteoarthritis (40). Osteoarthritis-induced chondrocyte apoptosis is prevented by microRNA-186 inhibition of the PI3K-AKT pathway via SPP1 (41). CXCL9 gene and protein expression were higher in proportion in the synovium of rheumatoid arthritis patients than that in osteoarthritis patients in a previous study (42). The same conclusions were drawn in a related study on the microarray analysis of gene expression in rheumatoid arthritis joints performed using cDNA microarrays and laser microdissection: CXCL9 and CXCL10 upregulation in the synovial lining associated with inflammation (43).

In GO and KEGG analyses, DEGs were enriched in immune response, chemokine-mediated signaling pathway, inflammatory response, chemokine activity, heparin-binding, antigen binding, cytokine-cytokine receptor interaction, chemokine signaling pathway, TNF signaling pathway, and NOD-like receptor signaling pathway. For receptor-chemokine interactions, the classical “two-site model” suggests that the globular core of chemokines serves as a docking structural domain of the receptor N-terminus (CRS1) that provides binding affinity, and chemokine N-terminus serves as a signaling trigger for the splice receptor binding pocket (CRS2) (44). One of the most important players in the disease process of two inflammatory joint diseases, rheumatoid arthritis and osteoarthritis, is the chemokine and chemokine receptor system (45). A recent study observed an increased expression of CC chemokine ligand 20 (a chemokine capable of binding to the CC chemokine receptor 6 expressed on Th17 cells) in inflamed joints of osteoarthritis patients (44). Cytokine-cytokine receptor interactions and rheumatoid arthritis are both triggered during the post-injury period (46). In addition, rheumatoid arthritis and osteoarthritis are caused by pro-inflammatory cytokines, which are released by the immune system in response to injury or inflammation (47). A local inflammatory process has been observed in osteoarthritic joints, and it is hypothesized that traumatic stimuli cause chondrocytes to produce cytokines and matrix metalloproteinases, which in turn lead to cartilage degradation in joints (48). A previous study confirmed the implication of inflammatory cytokines, TNF and IL-1, in the progression of rheumatoid arthritis and osteoarthritis (49).

Osteoarthritis is primarily a degenerative disease, rheumatoid arthritis is an autoimmune disease driven primarily by a significant inflammatory response involving the innate and adaptive immune systems (50). Inflammatory pathways of the innate immune system are activated in osteoarthritis because of the cellular stress and extracellular matrix degradation triggered by minor and major injury (51). A mosaic-like pattern of cytokine signaling and activation of molecular inflammatory pathways in the natural cells of intra-articular tissue pathogenesis contributes to the pathogenesis of osteoarthritis (52). Hence, The treatment of osteoarthritis through the modulation of the immune response is critical (53). In line with these findings, the GSEA analysis in this study revealed that osteoarthritis is predominantly linked to immune response activation. Furthermore, the GSEA data revealed a strong link to fatty acid metabolism. An important component of pathogenesis of osteoarthritis is the presence of intra-articular adipose tissue (54). Osteoarthritis is characterized by mitochondrial dysfunction and oxidative stress damage (55). For this reason, osteoblasts in human bone flap specimens from patients with osteoarthritis were studied and found to contain both in-situ active fatty acid oxidation and tricarboxylic acid metabolism (56). In osteoarthritis and other inflammatory disorders, fat pads contribute to the generation of adipokines and cytokines. This study and a previous study (57) demonstrated that adipokines regulate osteoarthritis and bone remodeling. This study’s GSEA analysis revealed that osteoarthritis was linked to graft-versus-host disease and Leishmania protozoa infection, which was one of the study’s findings.

The gene expression profile was subjected to Consensus Clustering by ConsensusClusterPlus package. In the omics analysis of large samples, it is often necessary to discuss the molecular typing of samples. The most common method in this paper is to cluster transcriptome, proteome and other data by Consensus Clustering. Finally, the samples can be divided into different clusters. There are obvious differences in transcriptome, proteome and other molecular patterns among the samples in different clusters, but the molecular patterns of the samples in each cluster are similar, thus realizing the purpose of molecular typing of large sample queues (58). For example, in the document “Proteogenomic landscape of squamous cell lung cancer”, based on the quantitative proteomic data of 108 lung squamous cell carcinoma samples, the author divided 108 tumor tissues into five molecular subtypes by consistent clustering, namely inflammatory subtype A, inflammatory subtype B, redox subtype A, redox subtype B and the molecular subtypes are obtained, the detailed subtype characteristics are described and discussed (59). A total of 494 lung squamous cell carcinoma samples were obtained from The Cancer Genome Atlas (TCGA) database. According to the proportion of immune cells identified in the samples, the tumor samples were analyzed by consistent clustering, and three subtypes were obtained. The subtypes I, II and III contained 17, 25 and 24 tumor samples, respectively. The clinical survival prognosis of subtype II is poor, while that of subtype I and III is good (60). In this study, we use consensus clustering to identify two immune patterns (clusterA and clusterB). ClusterA contains 34 samples and clusterB contains 35 samples. Through consistent cluster analysis and immune analysis, it is confirmed that there is significant heterogeneity between the two subgroups and the differential expression of related genes in the two subgroups. In the future, we will further study the specific mechanism of these genes in the immune microenvironment of osteoarthritis.

This study has some limitations. Firstly, we only studied the biological function of osteoarthritis in cells cultured in vitro, and further investigation of the unique processes is required. Second, although the aim of our study was to provide a landscape of mRNA and protein levels in osteoarthritis patients, we did not include more samples in our analysis (e.g., other rheumatic diseases). Future research should include a large number of samples and single cells. Additional histology (lipidomics, metabolomics, glycomics) and laboratory experiments may also improve the understanding of the conditions.

In summary, bioinformatics studies were performed to compare immune infiltration in the osteoarthritis and control groups. TCA1, TLR7, MMP9, CXCL10, CXCL13, HLA-DRA, and ADIPOQSPP1 were expressed at significantly higher levels in the IL-1β-induced group than in the osteoarthritis group, and our data also suggest that osteoarthritis may be associated with immune responses, chemokine-mediated signaling pathways, and inflammatory responses, and these studies improve our understanding of the development of osteoarthritis. The results of this study will aid in explaining the immune regulatory network of arthritis and inspire more effective therapeutic approaches. However, the functions of and associations between important genes and immune infiltration, as well as the significance of immune infiltration patterns in the development of osteoarthritis require further investigation.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/supplementary material.

Author Contributions

Data curation: KZ. Formal analysis: SN, YD. Methodology: SN, JQ, YD. Writing – original draft: XH. Writing – review and editing: XH, YD. All authors contributed to the article and approved the submitted version.

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.

Acknowledgments

The authors thank the staff of the Zhujiang Hospital, Southern Medical University.

Abbreviations

BP, Biological Process; MF, molecular function; CC, molecular function; GEO, Gene Expression Omnibus; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

References

1. Mi B, Wang J, Liu Y, Liu J, Hu L, Panayi AC, et al. Icariin Activates Autophagy via Down-Regulation of the NF-κb Signaling-Mediated Apoptosis in Chondrocytes. Front Pharmacol (2018) 9:605. doi: 10.3389/fphar.2018.00605

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Chang YH, Wu KC, Harn HJ, Lin SZ, Ding DC. Exosomes and Stem Cells in Degenerative Disease Diagnosis and Therapy. Cell Transplant (2018) 27:349–63. doi: 10.1177/0963689717723636

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Alahdal M, Duan L, Ouyang H, Wang D. The Role of Indoleamine 2,3 Dioxygenase 1 in Osteoarthritis. Am J Transl Res (2020) 12:2322–43.

PubMed Abstract | Google Scholar

4. Che X, Chen T, Wei L, Gu X, Gao Y, Liang S, et al. MicroRNA-1 Regulates the Development of Osteoarthritis in a Col2a1-Cre-ERT2/GFPfl/fl-RFP-miR-1 Mouse Model of Osteoarthritis Through the Downregulation of Indian Hedgehog Expression. Int J Mol Med (2020) 46:360–70. doi: 10.3892/ijmm.2020.4601

PubMed Abstract | CrossRef Full Text | Google Scholar

5. van der Kraan PM, Blaney Davidson EN, van den Berg WB. A Role for Age-Related Changes in TGFbeta Signaling in Aberrant Chondrocyte Differentiation and Osteoarthritis. Arthritis Res Ther (2010) 12:201. doi: 10.1186/ar2896

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Ripmeester EGJ, Timur UT, Caron MMJ, Welting TJM. Recent Insights Into the Contribution of the Changing Hypertrophic Chondrocyte Phenotype in the Development and Progression of Osteoarthritis. Front Bioeng Biotechnol (2018) 6:18. doi: 10.3389/fbioe.2018.00018

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Khosla S, Farr JN, Tchkonia T, Kirkland JL. The Role of Cellular Senescence in Aging and Endocrine Disease. Nat Rev Endocrinol (2020) 16:263–75. doi: 10.1038/s41574-020-0335-y

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Croft AP, Campos J, Jansen K, Turner JD, Marshall J, Attar M, et al. Distinct Fibroblast Subsets Drive Inflammation and Damage in Arthritis. Nature (2019) 570:246–51. doi: 10.1038/s41586-019-1263-7

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Tachmazidou I, Hatzikotoulas K, Southam L, Esparza-Gordillo J, Haberland V, Zheng J, et al. Identification of New Therapeutic Targets for Osteoarthritis Through Genome-Wide Analyses of UK Biobank Data. Nat Genet (2019) 51:230–6. doi: 10.1038/s41588-018-0327-1

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Wang J, Fan Q, Yu T, Zhang Y. Identifying the Hub Genes and Immune Cell Infiltration in Synovial Tissue Between Osteoarthritic and Rheumatoid Arthritic Patients by Bioinformatic Approach. Curr Pharm Des (2021) 28:497–509. doi: 10.2174/1381612827666211104154459

CrossRef Full Text | Google Scholar

11. Monteagudo S, Cornelis FMF, Aznar-Lopez C, Yibmantasiri P, Guns LA, Carmeliet P, et al. DOT1L Safeguards Cartilage Homeostasis and Protects Against Osteoarthritis. Nat Commun (2017) 8:15889. doi: 10.1038/ncomms15889

PubMed Abstract | CrossRef Full Text | Google Scholar

12. McAllister MJ, Chemaly M, Eakin AJ, Gibson DS, McGilligan VE. NLRP3 is a Potentially Novel Biomarker for the Management of Osteoarthritis. Osteoarthr Cartil (2018) 26:612–9. doi: 10.1016/j.joca.2018.02.901

CrossRef Full Text | Google Scholar

13. Xin X, Tan Q, Li F, Chen Z, Zhang K, Li F, et al. Potential Value of Matrix Metalloproteinase-13 as a Biomarker for Osteoarthritis. Front Surg (2021) 8:750047. doi: 10.3389/fsurg.2021.750047

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Orlov YL, Anashkina AA, Klimontov VV, Baranova AV. Medical Genetics, Genomics and Bioinformatics Aid in Understanding Molecular Mechanisms of Human Diseases. Int J Mol Sci (2021) 22(18):9962. doi: 10.3390/ijms22189962

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Wei C, Li J, Bumgarner RE. Sample Size for Detecting Differentially Expressed Genes in Microarray Experiments. BMC Genomics (2004) 5:87. doi: 10.1186/1471-2164-5-87

PubMed Abstract | CrossRef Full Text | Google Scholar

16. 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:R84. doi: 10.1186/ar4526

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Broeren MG, de Vries M, Bennink MB, Arntz OJ, Blom AB, Koenders MI, et al. Disease-Regulated Gene Therapy With Anti-Inflammatory Interleukin-10 Under the Control of the CXCL10 Promoter for the Treatment of Rheumatoid Arthritis. Hum Gene Ther (2016) 27:244–54. doi: 10.1089/hum.2015.127

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Broeren MG, de Vries M, Bennink MB, van Lent PL, van der Kraan PM, Koenders MI, et al. Functional Tissue Analysis Reveals Successful Cryopreservation of Human Osteoarthritic Synovium. PloS One (2016) 11:e0167076. doi: 10.1371/journal.pone.0167076

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Chakraborty S, Datta S, Datta S. Surrogate Variable Analysis Using Partial Least Squares (SVA-PLS) in Gene Expression Studies. Bioinformatics (2012) 28:799–806. doi: 10.1093/bioinformatics/bts022

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Bhattacharya S, Dunn P, Thomas CG, Smith B, Schaefer H, Chen J, et al. ImmPort, Toward Repurposing of Open Access Immunological Assay Data for Translational and Clinical Research. Sci Data (2018) 5:180015. doi: 10.1038/sdata.2018.15

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Stelzer G, Rosen N, Plaschkes I, Zimmerman S, Twik M, Fishilevich S, et al. The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses. Curr Protoc Bioinf (2016) 54:1.30.31–31.30.33. doi: 10.1002/cpbi.5

CrossRef Full Text | Google Scholar

22. 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:15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

23. 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:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Wilkerson MD, Hayes DN. ConsensusClusterPlus: A Class Discovery Tool With Confidence Assessments and Item Tracking. Bioinformatics (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Kanehisa M, Goto S. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res (2000) 28:27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Yu G, Wang LG, Han Y, He QY. Clusterprofiler: An R Package for Comparing Biological Themes Among Gene Clusters. Omics (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING V11proeinprotein Association Networks With Increased Coverage, Supporting Functional Discovery in Genome-Wide Experimental Datasets. Nucleic Acids Res (2019) 47:D607–13-d613. doi: 10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: A Softwareenvironment for Integrated Models of Biomolecular Interaction Networks. Genome Res (2003) 13:2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. ClueGO: A Cytoscape Plug-in to Decipher Functionally Grouped Gene Ontology and Pathway Annotation Networks. Bioinformatics (2009) 25:1091–3. doi: 10.1093/bioinformatics/btp101

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, et al. Systematic RNA Interference Reveals That Oncogenic KRAS-Driven Cancers Require TBK1. Nature (2009) 462:108–12. doi: 10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

31. 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:773–82. doi: 10.1038/s41587-019-0114-2

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Vincent TL. IL-1 in Osteoarthritis: Time for a Critical Review of the Literature. F1000Res (2019) 8:934. doi: 10.12688/f1000research.18831.1

CrossRef Full Text | Google Scholar

33. Xu Z, Ke T, Zhang Y, Guo L, Chen F, He W. Danshensu Inhibits the IL-1β-Induced Inflammatory Response in Chondrocytes and Osteoarthritis Possibly via Suppressing NF-κb Signaling Pathway. Mol Med (2021) 27(1):80. doi: 10.1186/s10020-021-00329-9

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Bowman S, Awad ME, Hamrick MW, Hunter M, Fulzele S. Recent Advances in Hyaluronic Acid Based Therapy for Osteoarthritis. Clin Transl Med (2018) 7:6. doi: 10.1186/s40169-017-0180-3

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Rosshirt N, Hagmann S, Tripel E, Gotterbarm T, Kirsch J, Zeifang F, et al. A Predominantthpolarization is Present in Synovial Fluid of End-Stage Osteoarthritic Knee Joints: Analysis of Peripheral Blood, Synovial Fluid and Synovial Membrane. Clin Exp Immunol (2019) 195:395–406. doi: 10.1111/cei.13230

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Jafri MA, Kalamegam G, Abbas M, Al-Kaff M, Ahmed F, Bakhashab S, et al. Deciphering the Association of Cytokines, Chemokines, and Growth Factors in Chondrogenic Differentiation of Human Bone Marrow Mesenchymal Stem Cells Using an Ex Vivo Osteochondral Culture System. Front Cell Dev Biol (2019) 7:380. doi: 10.3389/fcell.2019.00380

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Pandruvada SN, Yuvaraj S, Liu X, Sundaram K, Shanmugarajan S, Ries WL, et al. Role of CXCchemokine Ligand 13 in Oral Squamous Cell Carcinoma Associated Osteolysis in Athymic Mice. Int J Cancer (2010) 126:2319–29. doi: 10.1002/ijc.24920

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Kobayashi S, Watanabe T, Suzuki R, Furu M, Ito H, Ito J, et al. TGF-β Induces the Differentiation of Human CXCL13-Producing CD4(+) T Cells. Eur J Immunol (2016) 46:360–71. doi: 10.1002/eji.201546043

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Folkersen L, Brynedal B, Diaz-Gallo LM, Ramsköld D, Shchetynsky K, Westerlind H, et al. Integration of Known DNA, RNA and Protein Biomarkers Provides Prediction of Anti-TNF Response in Rheumatoid Arthritis: Results From the COMBINE Study. Mol Med (2016) 22:322–8. doi: 10.2119/molmed.2016.00078

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Liu R, Liu Q, Wang K, Dang X, Zhang F. Comparative Analysis of Gene Expression Profiles in Normal Hip Human Cartilage and Cartilage From Patients With Necrosis of the Femoral Head. Arthritis Res Ther (2016) 18:98. doi: 10.1186/s13075-016-0991-4

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Lin Z, Tian XY, Huang XX, He LL, Xu F. microRNA-186 Inhibition of PI3K-AKT Pathway via SPP1 Inhibits Chondrocyte Apoptosis in Mice With Osteoarthritis. J Cell Physiol (2019) 234:6042–53. doi: 10.1002/jcp.27225

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Chen HJ, Li Yim AYF, Griffith GR, de Jonge WJ, Mannens MMAM, Ferrero E, et al. Meta-Analysis of In Vitro-Differentiated Macrophages Identifies Transcriptomic Signatures That Classify Disease Macrophages In Vivo. Front Immunol (2019) 10:2887. doi: 10.3389/fimmu.2019.02887

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Yoshida S, Arakawa F, Higuchi F, Ishibashi Y, Goto M, Sugita Y, et al. Gene Expression Analysis of Rheumatoid Arthritis Synovial Lining Regions by cDNA Microarray Combined With Laser Microdissection: Up-Regulation of Inflammation-Associated STAT1, IRF1, CXCL9, CXCL10, and CCL5. Scand J Rheumatol (2012) 41:170–9. doi: 10.3109/03009742.2011.623137

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Zheng Y, Han GW, Abagyan R, Wu B, Stevens RC, Cherezov V, et al. Structure of CC Chemokine Receptor 5 With a Potent Chemokine Antagonist Reveals Mechanisms of Chemokine Recognition and Molecular Mimicry by HIV. Immunity (2017) 46:1005–1017.e5. doi: 10.1016/j.immuni.2017.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Hou SM, Chen PC, Lin CM, Fang ML, Chi MC, Liu JF. CXCL1 Contributes to IL-6 Expression in Osteoarthritis and Rheumatoid Arthritis Synovial Fibroblasts by CXCR2, C-Raf, MAPK, and AP-1 Pathway. Arthritis Res Ther (2020) 22:251. doi: 10.1186/s13075-020-02331-8

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Cheng XQ, Liang XZ, Wei S, Ding X, Han GH, Liu P, et al. Protein Microarray Analysis of Cytokine Expression Changes in Distal Stumps After Sciatic Nerve Transection. Neural Regener Res (2020) 15:503–11. doi: 10.4103/1673-5374.266062

CrossRef Full Text | Google Scholar

47. Schaible HG. Nociceptive Neurons Detect Cytokines in Arthritis. Arthritis Res Ther (2014) 16:470. doi: 10.1186/s13075-014-0470-8

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Kadam UT, Blagojevic M, Belcher J. Statin Use and Clinical Osteoarthritis in the General Population: A Longitudinal Study. J Gen Intern Med (2013) 28:943–9. doi: 10.1007/s11606-013-2382-8

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Stevens AL, Wheeler CA, Tannenbaum SR, Grodzinsky AJ. Nitric Oxide Enhances Aggrecan Degradation by Aggrecanase in Response to TNF-Alpha But Not IL-1beta Treatment at a Post-Transcriptional Level in Bovine Cartilage Explants. Osteoarthr Cartil (2008) 16:489–97. doi: 10.1016/j.joca.2007.07.015

CrossRef Full Text | Google Scholar

50. Paradowska-Gorycka A, Wajda A, Romanowska-Próchnicka K, Walczuk E, Kuca-Warnawin E, Kmiolek T, et al. Th17/Treg-Related Transcriptional Factor Expression and Cytokine Profile in Patients With Rheumatoid Arthritis. Front Immunol (2020) 11:572858. doi: 10.3389/fimmu.2020.572858

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Kraus VB, Blanco FJ, Englund M, Karsdal MA, Lohmander LS. Call for Standardized Definitions of Osteoarthritis and Risk Stratification for Clinical Trials and Clinical Use. Osteoarthr Cartil (2015) 23:1233–41. doi: 10.1016/j.joca.2015.03.036

CrossRef Full Text | Google Scholar

52. Cai A, Hutchison E, Hudson J, Kawashima Y, Komori N, Singh A, et al. Metabolic Enrichment of Omega-3 Polyunsaturated Fatty Acids Does Not Reduce the Onset of Idiopathic Knee Osteoarthritis in Mice. Osteoarthr Cartil (2014) 22:1301–9. doi: 10.1016/j.joca.2014.06.033

CrossRef Full Text | Google Scholar

53. Lin R, Deng C, Li X, Liu Y, Zhang M, Qin C, et al. Copper-Incorporated Bioactive Glass-Ceramics Inducing Anti-Inflammatory Phenotype and Regeneration of Cartilage/Bone Interface. Theranostics (2019) 9:6300–13. doi: 10.7150/thno.36120

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Eymard F, Pigenet A, Rose C, Bories A, Flouzat-Lachaniette CH, Berenbaum F, et al. Contribution of Adipocyte Precursors in the Phenotypic Specificity of Intra-Articular Adipose Tissues in Knee Osteoarthritis Patients. Arthritis Res Ther (2019) 21:252. doi: 10.1186/s13075-019-2058-9

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Chen P, Zheng L, Wang Y, Tao M, Xie Z, Xia C, et al. Desktop-Stereolithography 3D Printing of a Radially Oriented Extracellular Matrix/Mesenchymal Stem Cell Exosome Bioink for Osteochondral Defect Regeneration. Theranostics (2019) 9:2439–59. doi: 10.7150/thno.31017

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Karner CM, Long F. Glucose Metabolism in Bone. Bone (2018) 115:2–7. doi: 10.1016/j.bone.2017.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Peek V, Neumann E, Inoue T, Koenig S, Pflieger FJ, Gerstberger R, et al. Age-Dependent Changes of Adipokine and Cytokine Secretion From Rat Adipose Tissue by Endogenous and Exogenous Toll-Like Receptor Agonists. Front Immunol (2020) 11:1800. doi: 10.3389/fimmu.2020.01800

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Stewart PA, Welsh EA, Slebos RJC, Fang B, Izumi V, Chambers M, et al. Proteogenomic Landscape of Squamous Cell Lung Cancer. Nat Commun (2019) 10(1):3578. doi: 10.1038/s41467-019-11452-x

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Qiang H, Li J, Chang Q, Shen Y, Qian J, Chu T. Mining GEO and TCGA Database for Immune Microenvironment of Lung Squamous Cell Carcinoma Patients With or Without Chemotherapy. Front Oncol (2022) 12:835225. doi: 10.3389/fonc.2022.835225

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Li J, Xie L, Xie Y, Wang F. Bregmannian Consensus Clustering for Cancer Subtypes Analysis. Comput Methods Programs BioMed (2020) 189:105337. doi: 10.1016/j.cmpb.2020.105337

PubMed Abstract | CrossRef Full Text | Google Scholar