Sage Journals: Discover world-class research

Abstract

Objective

Although knee osteoarthritis (KOA) is a common disease, there is a lack of specific prevention and early treatment methods. Hence, this study aimed to examine the molecular changes occurring at different stages of KOA to elucidate the dynamic nature of the disease.

Design

Using a low-force compression model and analyzing RNA sequencing data, we identified molecular changes in the transcriptome of knee joint cartilage, including gene expression and molecular pathways, between the cellular changes and structural damage stages of KOA progression. In addition, we validated hub genes using an external dataset.

Results

Gene set enrichment analysis (GSEA) identified the following pathways to be associated with KOA: “B-cell receptor signaling pathway,” “cytokine-cytokine receptor interaction,” and “hematopoietic cell lineage.” Expression analysis revealed 585 differentially expressed genes, with 579 downregulated and 6 upregulated genes. Enrichment and clustering analyses revealed that the main molecular clusters were involved in cell cycle regulation and immune responses. Furthermore, the hub genes Csf1r, Cxcr4, Cxcl12, and Ptprc were related to immune responses.

Conclusions

Our study provides insights into the dynamic nature of early-stage KOA and offers valuable information to support the development of effective intervention strategies to prevent the irreversible damage associated with KOA, thereby addressing a major clinical challenge.

Keywords

animal model cartilage osteoarthritis RNA sequencing

Introduction

Osteoarthritis (OA) is a prevalent global condition characterized by joint cartilage degeneration and joint deformation.¹ Posttraumatic osteoarthritis (PTOA) is osteoarthritis that develops after joint injury. The knee is one of the most frequently injured joints in the body, possibly leading to knee osteoarthritis (KOA). Pain and deformity caused by KOA considerably impair lower limb function in patients. Conservative treatment and joint replacement surgery are the primary approaches to managing KOA. However, the prevention and treatment of KOA are limited to physical therapy and nonsteroidal anti-inflammatory drug (NSAID) use, with joint replacement surgery remaining the ultimate option for most patients.² Consequently, KOA is associated with a high disability rate, which severely decreases the quality of life of patients and imposes a substantial economic burden on families and society.^3,4 To develop effective strategies for KOA prevention and treatment, it is necessary to investigate the underlying mechanisms of its early development.

The onset of KOA pathology may be slow,⁵ with the earliest stages limited to only cellular changes (de Windt et al., 2013)⁶—“pre-osteoarthritis” (pre-OA) stage⁷—making it crucial to employ molecular biology and bioinformatic techniques to elucidate these molecular changes. Animal models are an effective approach for studying KOA progression. However, most studies have focused on inducing rapid knee joint destruction, reflecting end-stage OA tissue,^8-11 with limited attention paid to the early progression of KOA before the occurrence of irreversible structural damage. We previously developed a rat model in which KOA was induced through a single session of low-force cyclic compression.¹² Our observations revealed that cysts and cavities in the tibial cartilage appeared 6 months after compression, while fissures in the cartilage were observed after 1 year of compression. We believe that this model can provide valuable insights into the progression of pre-structural damage in pre-OA patients. Therefore, this study aimed to explore the molecular changes at the transcriptomic level, from cellular changes to irreversible structural damage, including changes in gene expression, molecular pathways, and interaction networks.

Methods

Animals and Treatments

Healthy male Wistar rats at 12 weeks of age (provided by SHIMIZU Laboratory Supplies Co., Ltd., Kyoto, Japan) were used. The rats were housed in a standard environment at a temperature of 24°C, humidity of 55% ± 5%, good ventilation, and ample space for activity. The rats were provided food (standard solid pellet feed) and water. The rats in the compression group were subjected to low-force compression as reported previously.^13,14 Briefly, a single session of cyclic compression, which did not cause structural damage to the right knee joint, was applied. This was followed by normal feeding (adequate standard solid pellet feed, water, and exercise space) until the desired time points (8 rats allocated for each time point). The rats in the control group underwent identical anesthesia procedures but did not undergo compression. They were collected at the same time points as in the compression group (8 rats allocated for each time point). Previously,¹² we defined the pre-structural damage stage as 12 weeks post-compression (+12 weeks group) and the structural damage stage as 6 months post-compression (+6 months group).

Histology

The rats were euthanized by overdose of isoflurane inhalation followed by exsanguination. Tissue blocks were collected from the right knee joint of each rat. For histological observations, tissue blocks were fixed in 4% paraformaldehyde, decalcified in 10% ethylenediaminetetraacetic acid, and embedded in paraffin. Joint tissue blocks were sectioned in the sagittal plane. Sections were mounted on glass slides and deparaffinized. For Safranin O/Fast Green staining, the sections were stained with 0.02% Fast Green followed by 0.1% Safranin O. Tibial cartilage damage was quantified using a modified Mankin scoring system.¹⁵

RNA-Sequencing and Bioinformatic Analysis

Cartilage samples were collected from the tibial plateau of the right knee, immediately frozen in liquid nitrogen, and stored at −80°C. Total RNA was extracted using an RNeasy Plus Universal Mini Kit (Qiagen, Valencia, CA, USA) following the manufacturer’s protocol. Each RNA sample for total RNA sequencing (RNA-Seq) was extracted from 2 individuals treated under the same conditions simultaneously and pooled to ensure sufficient RNA quantity while reducing the impact of individual variability. Analyses were performed in 4 replicates (N = 4 pooled RNA samples per group and time point). RNA quality assessment, cDNA library construction, and RNA sequencing were conducted at SignAC (Single-cell Genome Information Analysis Core at WPI-ASHBi, Kyoto University, Kyoto, Japan). Samples with an RNA integrity number (RIN) greater than 7, indicating intact RNA, were used for RNA-Seq. Sequencing was performed using NovaSeq6000 with 151-bp paired-end reads.

After sequencing, the SnakePipes (v. 2.7.0) mRNA-Seq pipeline was used for data processing and analysis. Low-quality reads, adapter contamination, and high-content unknown nucleotide (N) noise reads were filtered out using Fastp (v 0.23.2). The remaining “clean reads” were aligned to the reference genome (NCBI mRatBN7.2) using STAR (v 2.7.10a). Low expression data (molecules with 0 values greater than 50% across all samples) and genes that had not yet been specifically named or annotated were further filtered out. Sequence data from each group were analyzed. Principal component analysis (PCA) (version 3.6.0), k-means cluster analysis, and gene set enrichment analysis (GSEA) (version 4.8.1) were performed using the transformed data. DEseq2¹⁶ was used to detect differentially expressed genes (DEGs). DEG analysis was conducted between the +6 months and +12 weeks groups using a threshold of |logFC| > 1 and false discovery rate (FDR) < 0.05, for both the compression and control groups. Functional annotation was performed using the ClusterProfiler software (version 4.8.1). Protein-protein interactions (PPIs) were analyzed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING)¹⁷ and Cytoscape (v. 3.9.1, Institute of Systems Biology, Seattle, WA).¹⁸

Validation of the Hub Genes

External validation was conducted to evaluate the capacity of the hub genes in distinguishing between cellular changes and structural damage stage using data from Gene Expression Omnibus (GEO) databases (GSE28958, GSE42295, and GSE103416). Table 1 lists the characteristics of the datasets. To create an integrated GEO dataset, batch effects were removed,^19,20 and the samples were grouped based on known histological information.^21,22 The integrated dataset was then subjected to receiver-operating characteristic (ROC) analysis using pROC (v. 1.17.0.1) to obtain the area under the curve (AUC).

Table 1.

Descriptive Statistics.

Data Number	Platform	KOA Model	Early Injury Group	Damage Group	Species
GSE28958	GPL6247	Monosodium iodoacetate	OA + 5 days, N = 3	OA + 21 days, N = 3	Rattus norvegicus
GSE42295	GPL1355	Surgical destabilization	OA + 2 weeks, N = 3	OA + 8 weeks, N = 3	Rattus norvegicus
GSE103416	GPL17117	Monosodium iodoacetate	OA + 2 days, N = 4	OA + 28 days, N = 4	Rattus norvegicus

KOA = knee osteoarthritis; OA = osteoarthritis.

Overall Workflow

The overall workflow involved the RNA-Seq of rat cartilage samples at different time points after being subjected to low-force cyclic compression. The 12-week time point (+12 weeks), characterized by minimal observable structural cartilage damage and early injury manifesting as cellular changes, was designated as the early injury stage. Whereas the 6-month time point (+6 months), characterized by the presence of intra-articular cysts and cavity damage, was designated as the damage stage. Differential analysis was conducted between the compressed cartilage samples at +12 weeks post-compression and those at +6 months post-compression to elucidate the progression of cartilage injuries following compression. Similarly, cartilage samples from the control group underwent differential analysis at the same 2 time points to elucidate the intrinsic background changes in the cartilage over time. The data analysis process included quality control, genome alignment, gene expression quantification, differential analysis, clustering analysis, GSEA, gene function enrichment, PPI network analysis, and external dataset validation ( Fig. 1 ).

Figure 1.

Study flow chart. PCA = principal component analysis; GSEA = gene set enrichment analysis; OA = osteoarthritis; DEGs = differentially expressed genes; KEGG = Kyoto Encyclopedia of Genes and Genomes; GO = Gene Ontology; ROC = receiver-operating characteristic; PPI = protein-protein interaction.

Statistical Analysis

Tissue histology data were statistically analyzed using the GraphPad software (v. 9.4.1). The non-parametric Mann-Whitney U test was employed to compare differences among groups, and differences were considered statistically significant at P < 0.05.

Results

Morphology of Cartilage at Different Time Points

Figure 2A displays the macroscopic tissue images. The cartilage surface of the +12 weeks group exhibited slight granularity compared with the uncompressed cartilage, whereas visible cartilage surface ulcers were observed in the +6 months group. Figure 2B shows the histological observations. In the +12 weeks group rats, the cartilage showed only features of reduced staining and hypertrophic nuclei, whereas in the +6 months group rats, the cartilage exhibited evident intra-articular cysts and cavities in the middle layer. Histological scoring indicated significant KOA progression at +6 months compared to that at +12 weeks ( Fig. 2C ).

Figure 2.

Morphological observation and analysis. (A) Macroscopic observation of the tibial cartilage. “M” represents the medial side of the tibia and “L” represents the lateral side. The red arrows indicate cartilage surface ulcers. (B) Histological observation of the tibial cartilage stained with Safranin O and Fast Green. The yellow arrows indicate structural damage to the cartilage. (C) Results of modified Mankin scoring. *P < 0.05.

Gene Expression Profiles Among the Groups

A total of 33,294 genes were identified across all samples and subjected to further analysis after screening. The median, interquartile range, and maximum and minimum values of gene expression were comparable among samples ( Fig. 3A ). Since all samples were derived from OA cartilage, a relatively high degree of correlation among all samples was expected among them. Nevertheless, correlation analysis revealed stronger correlations within the +6 months group samples ( Fig. 3B ). PCA demonstrated a distinction between the sample clusters at the 2 time points, indicating differential gene expression between the groups ( Fig. 3C ). GSEA provided enrichment results for the entire expressed gene set between the 2 time points. The GSEA revealed enrichment related to well-established immune processes such as the “B-cell receptor signaling pathway,” “cytokine-cytokine receptor interaction,” and “hematopoietic cell lineage.”

Figure 3.

Analysis of gene expression profiles of the +12 weeks and +6 months groups. (A) Boxplots of gene probe expression. (B) Heatmap of correlation among the samples. Numerical representation of the degree of correlation between samples. (C) Principal component analysis. (D) Gene set enrichment analysis. The figure indicates the top 3 enriched processes. OA = osteoarthritis.

K-means analysis was performed on the 2,000 most variably expressed genes ( Fig. 4A ). The number of clusters was set to 4, based on the elbow method, with cluster C containing the highest number of genes. Further enrichment analyses were conducted for each cluster ( Fig. 4B ). The enrichment analysis results were similar to the GSEA results, with the expression of most genes being associated with immune-related processes.

Figure 4.

K-means cluster analysis and enrichment results. (A) A heatmap of the 2,000 most variably expressed genes. (B) Enrichment of clusters. Each color represents a distinct cluster. The size of the dot indicates the number of enriched genes. The connections represent clustering relationships. Each entry comprises a P-value and term name.

Identification of Differentially Expressed Genes

In total, 585 DEGs were identified between the +6 months and +12 weeks groups (Supplement 1). Among these, 579 genes were downregulated and 6 genes were upregulated ( Fig. 5A ). The expression of the widely known OA-associated markers Adamts9 and matrix metalloproteinase 13 (MMP13), and some other MMP family members, col1a1, and col1a2 were significantly downregulated in the +6 months group compared to that in the +12 weeks group. Meanwhile, col2a1 expression was not significantly different between the groups. A heat map of the top 50 most significant DEGs was generated to visualize the expression patterns of the genes in the cartilage samples at the 2 time points ( Fig. 5B ). Lilrb3a, Anpep, Cd44, Igfbp3, and Ctsk were the top 5 downregulated genes, with logFC magnitudes ranging from −2.14 to −1.13. Kcnj4, Csmd2, Kcnk3, Amigo2, Grtp1, and Serpina12 were the 6 upregulated genes, with logFC magnitudes ranging from 1.01 to 3.42. Downregulated DEGs primarily participated in “immune system processes” and “immune response,” and upregulated DEGs were involved in “potassium ion transport” ( Fig. 5C ).

Figure 5.

DEG analysis. (A) Volcano plot of DEGs. The green dots represent downregulated genes and the red dots represent upregulated genes. (B) Heatmap of the top 50 most significant DEGs. The color key represents the Z-score. (C) DEG pathway network. Two pathways (nodes) are connected if they share 30% or more genes. Green and red represent downregulated and upregulated pathways. Darker nodes represent more significantly enriched gene sets. Bigger nodes represent larger gene sets. Thicker edges represent more overlapped genes. DEG = differentially expressed gene; OA = osteoarthritis.

Enrichment of DEGs

All DEGs were annotated using Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. GO enrichment focused on gene functions and processes. In the GO term enrichment analysis, the top 3 significantly enriched terms in the categories of biological process (BP), cellular component (CC), and molecular function (MF) were “leukocyte migration,” “myeloid leukocyte activation,” and “leukocyte proliferation”; “external side of plasma membrane,” “integrin complex,” and “protein complex involved in cell adhesion”; and “cell adhesion molecule binding,” “cytokine binding,” and “integrin binding,” respectively ( Fig. 6A-C ). KEGG enrichment focuses on pathways and functional modules. In the KEGG pathway enrichment analysis, the top 3 significantly enriched pathways were “hematopoietic cell lineage,” “platelet activation,” and “B-cell receptor signaling pathway” ( Fig. 6D ). The GO and KEGG enrichment analysis results were immune-related.

Figure 6.

DEG enrichment. (A-C) GO enrichment analysis in the categories of BP, CC, and MF. (D) KEGG pathway enrichment analysis. The size of the dots indicates the number of genes involved and the color key indicates the −log10 transformed P-values. The terms are sorted in a descending order based on the gene ratio. DEG = differentially expressed gene; BP = biological process; CC = cellular component; MF = molecular function; GO = Gene Ontology; KEGG = Kyoto Encyclopedia of Genes and Genomes.

Removal of the Background Variation

To eliminate the influence of background changes during cartilage tissue development over time, we performed RNA-Seq and data acquisition on normal control cartilage at the same time points to identify control-DEGs. We identified 62 downregulated and 40 upregulated control-DEGs between the control +6 months and control +12 weeks ( Fig. 7A ). Among these, 39 were downregulated ( Fig. 7B ), and 1 was upregulated as background-related DEGs ( Fig. 7C ). These background-related DEGs were found to be shared with OA-DEGs, indicating that, over time, normal cartilage undergoes molecular expression changes. We excluded these 40 background-related DEGs from OA-DEGs, which resulted in a set of DEGs referred to as “de-background DEGs” (Supplement 1).

Figure 7.

Identification of background variation. (A) Volcano plot of control DEGs. The green dots represent downregulated genes and the red dots represent upregulated genes. (B-C) The Venn diagram of OA-DEGs and control-DEGs. The overlapping region contains background DEGs. DEGs = differentially expressed genes; OA = osteoarthritis.

PPI Network and Identification of Cluster

We constructed a PPI network to explore the interactions between the proteins encoded by the DEGs. Subsequently, we used a Markov clustering algorithm to cluster the molecules, resulting in 114 clusters (Supplement 2). Seven of these clusters contained 10 or more genes ( Fig. 8A ). Among these clusters, Cluster 1 and Cluster 2, with 77 and 54 members, respectively, aggregated in the network. The internal interactions within Clusters 1 and 2 were visualized using Cytoscape. In Cluster 1, the number of nodes was 77, the number of edges was 1,698, and the average node degree was 44.1 ( Fig. 8B ). In Cluster 2, the number of nodes was 54, the number of edges was 342, and the average node degree was 12.7 ( Fig. 8C ).

Figure 8.

Construction of a PPI network and generation of Markov clustering. (A) PPI network of de-background differentially expressed genes. Clusters including 10 or more members are shown. Nodes of the same color represent proteins belonging to the same cluster. The edges in the network represent protein-protein associations and the width of the edges reflects the confidence level of the associations. (B-C) PPI network of Clusters 1 and 2. Nodes are represented by circles, and the size and color intensity of the circles indicate the degree of interaction. PPI = protein-protein interaction.

Cluster-DEG Functional Annotation

Enrichment analyses of Clusters 1 and 2 were performed, annotating them to the GO and KEGG databases. In the GO term enrichment analysis, Cluster 1 genes were found to be associated with cell cycle, with the most significantly enriched BP, CC, and MF terms being “chromosome segregation,” “spindle,” and “microtubule binding,” respectively ( Fig. 9A-C ). Cluster 2 genes showed significant enrichment in “leukocyte migration,” “external side of the plasma membrane,” and “immune receptor activity” ( Fig. 10A-C ). The top enriched KEGG pathway of Cluster 1 genes was “cell cycle” ( Fig. 9D ), whereas that of Cluster 2 genes was “hematopoietic cell lineage” ( Fig. 10D ), suggesting distinct molecular patterns between the clusters. These findings are consistent with the results of GSEA and K-means clustering, confirming the immune-related nature of Cluster 2.

Figure 9.

Enrichment analysis of cluster 1 DEGs. (A-C) BPs, CCs, and MFs enriched by DEGs. (D) KEGG pathway enrichment. Different colors denote different terms, and connecting lines represent the participating molecules. DEGs = differentially expressed genes; BP = biological process; CC = cellular component; MF = molecular function; KEGG = Kyoto Encyclopedia of Genes and Genomes.

Figure 10.

Enrichment analysis of cluster 2 DEGs. (A-C) BPs, CCs, and MFs enriched by DEGs. (D) KEGG pathway enrichment. Different colors indicate different terms, and connecting lines represent the participating molecules. DEGs = differentially expressed genes; BP = biological process; CC = cellular component; MF = molecular function; KEGG = Kyoto Encyclopedia of Genes and Genomes.

Identification of Hub Genes

Since many analytical results point toward immune-related processes, we aimed to further identify hub immune-related genes. Functional module construction using the MCODE plugin in Cytoscape was performed using Cluster 2 DEGs. The results revealed the presence of 3 modules with a top module score of 11.529, comprising 18 nodes and 98 edges. Subsequently, the top 10 hub genes were selected using the 12 different topological methods provided by the CytoHubba plugin in Cytoscape. Among these modules, 10 genes (Ptprc, Tlr7, Csf1r, Cd4, Itgam, Sell, Cd44, Cxcr4, Itgax, and Cxcl12) appeared in over half of the modules ( Table 2 ).

Table 2.

Twelve Topologically Selected Genes.

MCC	DMNC	MNC	Degree	EPC	Bottleneck	Eccentricity	Closeness	Radiality	Betweenness	Stress	Clustering Coefficient
Ptprc	Cyth4	Ptprc	Cd4	Ptprc	Ptprc	Itgam	Cd4	Cd4	Cd4	Cd4	Cnr2
Tlr7	Rgs18	Cd4	Ptprc	Cd4	Cd4	Ptprc	Ptprc	Ptprc	Ptprc	Ptprc	Atp8a1
Fyb	Fgl2	Itgam	Itgam	Itgam	Sell	Cxcl12	Itgam	Itgam	Itgam	Itgam	Defa5
Rgs18	Fyb	Csf1r	Csf1r	Csf1r	Cxcl12	Cd38	Csf1r	Csf1r	Csf1r	Csf1r	Np4
Plek	Plek	Sell	Sell	Sell	Itgam	Gfap	Sell	Sell	Prom1	Cxcr4	Pag1
Cyth4	Prom1	Tlr7	Tlr7	Tlr7	Csf1r	Cxcr4	Tlr7	Cxcr4	Cxcr4	Tlr7	Dntt
Csf2rb	Ikzf1	Cxcr4	Cxcr4	Cxcr4	Tlr7	Cd44	Cxcr4	Tlr7	Sell	Sell	Fgl2
Csf1r	Spn	Itgax	Cd44	Itgax	Siglec5	Cd4	Cd44	Cd44	Cd44	Cd44	Tlr8
Ikzf1	Wdfy4	Cd44	Itgax	Cd68	Wdfy4	Sell	Itgax	Itgax	Tlr7	Prom1	Ccl6
Wdfy4	Pld4	Cxcl12	Cxcl12	Cd44	Prom1	Mpo	Cxcl12	Cxcl12	Itgax	Cd68	Cyth4

MCC = Maximal Clique Centrality; DMNC = Density of Maximum Neighborhood Component; MNC = Maximum Neighborhood Component; EPC = Edge Percolated Component.

Validation of the Hub Genes

The ability of the hub genes to discriminate between early and late stages of OA in previously published datasets was evaluated. Three GEO datasets were retrieved and integrated. After batch effect removal, the median expression levels were consistent across datasets ( Fig. 11A ). ROC analysis was performed to validate hub genes between the 2 time points. Excluding Tlr7 and Itgax, which were not included in the integrated dataset, the other genes showed good classification performance ( Fig. 11B ). The AUC values for Csf1r, Cxcr4, Cxcl12, and Ptprc were above 0.8, and the 95% confidence intervals were higher than 0.6.

Figure 11.

Validation of hub genes. (A) Integration results of 3 datasets. Three different colors represent different datasets. On the right, the expression levels of genes before and after dataset integration are shown. (B) Receiver-operating characteristic curves for 8 characterized genes to independently distinguish between the early injury and damage stages. CI = confidence interval.

Discussion

This study reports the molecular changes and characteristics of the rat knee joint cartilage transitioning from cellular to structural OA damage. Throughout this progression, most DEGs were downregulated, with immunological processes such as the “B-cell receptor signaling pathway” prominently featured. Among these DEGs, a clear cluster emerged: one cluster participated in a series of cell cycle–related BPs, while the other engaged in immune-related processes exemplified by “leukocyte migration.” Importantly, these molecules reflected a decrease in their gene expression compared to early disease, surpassing the levels that could be attributed solely to aging. This suggests a progression of pathology from cellular to structural alterations.

The study of OA progression within a time series is crucial because failure to grasp the progression of the disease can have implications for disease management. However, research on this topic is limited. Most studies have focused primarily on comparing end-stage OA tissues with their normal counterparts. Previously reported upregulated genes in late-stage OA tissues, including cartilage, synovium, and subchondral bone, have exhibited enrichment in pathways and processes associated with protein hydrolysis, extracellular matrix degradation, and collagen breakdown metabolism, whereas downregulated genes are related to cellular proliferation and the responsiveness of cells to stimuli.^9,11,23,24 Our results indicate that this downregulation was also significant in OA cartilage over time as the disease progressed from cellular changes to structural damage, suggesting that OA has a relatively long disease course and exhibits distinct BPs at various developmental stages. Early-stage disease is more likely to manifest as initial cellular changes in response to stimuli, whereas late-stage disease is dominated by tissue remodeling characterized by structural alterations in the joint.^25,26 In future studies on OA treatment approaches, it may be important to consider the effects of the dynamic nature of the disease itself. Neglecting this aspect may result in incorrect identification or overestimation of treatment effects.

Chondrocytes are the only cell type found in healthy cartilage and are responsible for producing collagen and proteoglycans to maintain the cartilage matrix.²⁷ In normal joint cartilage, chondrocytes remain quiescent, have ceased differentiation, and maintain an unmineralized matrix.^28,29 However, during OA, chondrocytes frequently undergo hypertrophy, followed by matrix alterations leading to cartilage loss. Therefore, the pathogenesis of OA is thought to recapitulate developmental processes.³⁰ In our study, the KEGG pathway and GO term enrichment analyses revealed numerous important terms associated with cell cycle and mitosis-related processes in Cluster 1 of DEGs. Within Cluster 1, genes involved in the cell cycle and spindle formation, such as Cdk1, Ccnb1, Ccnb2, Aurkb, and Bub1b, were identified. The transcriptional products of these genes exhibit a high degree of interaction with other proteins in the cluster. In addition, Ccnb1, Aurkb, and Bub1b in Cluster1 are involved in chromosome segregation. Cdk1 is a crucial cell mitosis regulator that regulates the transition from the G2 phase to the M phase of the cell cycle by interacting with the cell cycle proteins cyclin A1 and cyclin B1.^31,32 Cdk1, Ccnb1, and Ccnb2 of Cluster 1 are involved in the p53 signaling pathway and cellular senescence. Alterations in the expression of these proteins profoundly affect the G2 cell cycle arrest. When their expression is decreased, the cell cycle is halted in the G2/M phase, leading to the inability of cells to undergo mitosis and cell death.³¹ p53 plays a critical role in regulating apoptosis^33,34 and autophagy,³⁵ which are typical pathological changes in OA.³⁶ A previous study has indicated that p53 deletion rescues chondrocyte apoptosis, but fails to enhance chondrocyte proliferation.³⁷ Aurkb is involved in regulating chromosome alignment and separation during mitosis and meiosis by binding to microtubules. Bub1b is expected to play a role in the mid-to-late transition of the mitotic cell cycle and protein localization to chromosomes and kinetochore regions. Aurkb and Bub1b have been identified as potential target genes in microarray-based studies on OA.^38,39

In our study, we observed significant involvement of immune-related mechanisms during the transition from the stage with only cellular changes to the stage with the appearance of structural cystic lesions in the cartilage. In this study, GSEA of all genes expressed at both injury time points demonstrated enrichment of the B-cell receptor signaling pathway. B cells originate from pluripotent hematopoietic stem cells and are part of the adaptive immune system. B-cell response dysregulation can ultimately lead to autoimmunity,⁴⁰ making B cells a common target for treating inflammatory arthritis and systemic autoimmune diseases. Recent studies have shown changes in B lymphocyte responses in patients with OA. Chronic B-cell activation occurs in OA, and these changes cannot be solely explained by age.⁴¹ Consistent with our GSEA results, enrichment analyses of k-means clustering and cluster 2 DEGs indicated that immune response–related processes exhibit a significant pattern of transition from the early injury stage to the damage stage. The Csf1r-Cxcl12/Cxcr4 axis found in Cluster 2 is involved in leukocyte-related processes such as leukocyte migration and leukocyte chemotaxis. Csf1r promotes macrophage differentiation, activation, and osteoclastogenesis.⁴² Csf1r has been reported to participate in various processes, including osteoclast differentiation regulation and bone resorption modulation.⁴³ Studies have shown that inhibiting Csf1r expression has a protective effect on bone and cartilage in experimental arthritis.^44,45 The Csf1r-Cxcl12/Cxcr4 axis may play a crucial role in the inflammatory response of OA, and our validation using external datasets also demonstrated that these molecules exhibit the best classification performance. Ptprc is an important antigen receptor signaling regulator in T cells and B cells.⁴⁶ It plays a role in modulating the JAK/STAT pathway targeted in inflammatory diseases.⁴⁷ A previous study has indicated that the loss of Ptprc leads to an increase in JAK/STAT signaling.⁴⁸

The distinctive feature of the compression model utilized in this study lies in its relatively mild dose, which does not immediately induce substantial structural changes in the joint, but rather triggers only minor cellular alterations. This sets it apart from many other commonly employed PTOA models that induce rapid OA onset and progression through high-energy physical impacts or surgical procedures. For example, rats underwent anterior cruciate ligament transection and partial medial meniscectomy, followed by three 30-minute sessions of forced exercise each week. After only 4 weeks, disordered cartilage fiber alignment became evident.²¹ Macroscopic cartilage damage was observed, and a significant reduction in cartilage thickness was noted within 9 days of injecting monoiodoacetic acid into the rat knee joint. After 21 days, nearly the entire cartilage layer had disappeared.²² These rapid models help shorten the experimental timeline but cannot accurately reflect OA. A slower disease progression may provide a favorable window of opportunity for research into the prevention and intervention of PTOA. It is worth noting that despite the use of various methods to induce the rapid onset and progression of OA, the hub genes we selected performed well in distinguishing between the early injury and damage stages of KOA. This corroborates the existence of a stage, regardless of the triggering event, where cartilage function remains intact; however, BPs at the cellular level have already commenced joint degradation.⁷ This stage represents the optimal opportunity for implementing intervention measures to prevent irreversible structural damage.

This study had several limitations. First, it was conducted with only male animals, and further investigations are warranted to assess the impact of sex and validate the findings in human samples. Second, owing to experimental constraints, we compared different individuals at different time points of development, and additional experimental methods should be designed to include self-control within the same individual. Finally, although our study aimed to provide an overview of the changes in mRNA and protein levels during the transition from the early injury to the damage stage, we did not include numerous samples in our analysis. Furthermore, since all samples are from the early injury to the damage stage and given that our sample size was small, individual variation between samples likely influenced the correlation analysis and PCA results. Our DEG threshold was set conservatively to avoid false negatives. Future studies should use larger sample sizes. Additional omics studies, single-cell analyses, and laboratory experiments can improve our understanding of this disease.

In summary, observing and comparing the different stages of KOA according to the duration of OA is necessary. In our study, we observed downregulation of genes from the early injury stage to the damage stage of KOA. The molecules encoded by these genes exhibit 2 distinct patterns: one is involved in the transition of chondrocytes from a quiescent state to a differentiated state, accompanied by cell death, and the other is involved in immune-related processes, indicating immune dysregulation within the cartilage microenvironment. Csf1r, Cxcr4, Cxcl12, Ptprc, Sell, Itgam, Cd44, and Cd4 may be key candidate genes for immune-related patterns. However, further research is required to elucidate the functional roles of these important genes and their causal relationships with immune responses. These results will be useful for future research on the prevention and treatment of KOA.

Supplemental Material

sj-xlsx-1-car-10.1177_19476035231213174 – Supplemental material for Analysis of Molecular Changes and Features in Rat Knee Osteoarthritis Cartilage: Progress From Cellular Changes to Structural Damage

Supplemental material, sj-xlsx-1-car-10.1177_19476035231213174 for Analysis of Molecular Changes and Features in Rat Knee Osteoarthritis Cartilage: Progress From Cellular Changes to Structural Damage by Zixi Zhao, Akira Ito, Hiroshi Kuroki and Tomoki Aoyama in CARTILAGE

Supplemental Material

sj-xlsx-2-car-10.1177_19476035231213174 – Supplemental material for Analysis of Molecular Changes and Features in Rat Knee Osteoarthritis Cartilage: Progress From Cellular Changes to Structural Damage

Supplemental material, sj-xlsx-2-car-10.1177_19476035231213174 for Analysis of Molecular Changes and Features in Rat Knee Osteoarthritis Cartilage: Progress From Cellular Changes to Structural Damage by Zixi Zhao, Akira Ito, Hiroshi Kuroki and Tomoki Aoyama in CARTILAGE

Footnotes

Author Contributions

All authors were involved in drafting and revising the manuscript and approved the final version for publication. Study conception and design: Zixi Zhao, Akira Ito, Tomoki Aoyama, and Hiroshi Kuroki. Analysis and interpretation of data: Zixi Zhao and Akira Ito. Zixi Zhao drafted the manuscript.

Acknowledgments and Funding

We thank the Single-cell Genome Information Analysis Core (SignAC) at WPI-ASHBi, Kyoto University, for their support. The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by JST SPRING, grant JPMJSP2110 (ZZ), JSPS KAKENHI, grants JP21H03302 (AI), JP21K19709 (HK), JP23H03245 (HK), and JP18H03129 (HK).

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Ethical Approval

This study was approved by the Animal Research Committee of Kyoto University (approval number: Medkyo21082).

ORCID iDs

Zixi Zhao

Akira Ito

Data Availability

Computational methods were implemented in the R package, which can be downloaded from various websites. All external dataset data can be found in GEO under the accession code. All other data are available in the manuscript and supplementary files.

Supplementary material for this article is available on the Cartilage website at .

References

Martel-Pelletier

Barr

Cicuttini

Conaghan

Cooper

Goldring

, et al. Osteoarthritis. Nat Rev Dis Primers. 2016;2:16072. doi:10.1038/nrdp.2016.72.

Evans

Walker

Evans

Blom

Sayers

Whitehouse

. How long does a knee replacement last? A systematic review and meta-analysis of case series and national registry reports with more than 15 years of follow-up. Lancet. 2019;393:655-63. doi:10.1016/S0140-6736(18)32531-5.

Peat

Thomas

. Osteoarthritis year in review 2020: epidemiology & therapy. Osteoarthritis Cartilage. 2021;29:180-9. doi:10.1016/j.joca.2020.10.007.

Safiri

Kolahi

Smith

Hill

Bettampadi

Mansournia

, et al. Global, regional and national burden of osteoarthritis 1990-2017: a systematic analysis of the Global Burden of Disease Study 2017. Ann Rheum Dis. 2020;79(6):819-28. doi:10.1136/annrheumdis-2019-216515.

van Saase

van Romunde

Cats

Vandenbroucke

Valkenburg

. Epidemiology of osteoarthritis: Zoetermeer survey. Comparison of radiological osteoarthritis in a Dutch population with that in 10 other populations. Ann Rheum Dis. 1989;48:271-80. doi:10.1136/ard.48.4.271.

de Windt

Vonk

Brittberg

Saris

DBF

. Treatment and prevention of (early) osteoarthritis using articular cartilage repair-fact or fiction? A systematic review. Cartilage. 2013;4:5S-12S. doi:10.1177/1947603513486560.

Ryd

Brittberg

Eriksson

Jurvelin

Lindahl

Marlovits

, et al. Pre-osteoarthritis. Cartilage. 2015;6:156-65. doi:10.1177/1947603515586048.

Aigner

Fundel

Saas

Gebhard

Haag

Weiss

, et al. Large-scale gene expression profiling reveals major pathogenetic pathways of cartilage degeneration in osteoarthritis. Arthritis Rheum. 2006;54:3533-44. doi:10.1002/art.22174.

Davidson

Waters

Kevorkian

Darrah

Cooper

. Expression profiling of metalloproteinases and their inhibitors in synovium and cartilage. Arthritis Res Ther. 2006;8:R124. doi:10.1186/ar2013.

10.

Karlsson

Dehne

Lindahl

Brittberg

Pruss

Sittinger

, et al. Genome-wide expression profiling reveals new candidate genes associated with osteoarthritis. Osteoarthritis Cartilage. 2010;18:581-92. doi:10.1016/j.joca.2009.12.002.

11.

Chou

Lee

Song

Chuang

Kuo

, et al. Direct assessment of articular cartilage and underlying subchondral bone reveals a progressive gene expression change in human osteoarthritic knees. Osteoarthritis Cartilage. 2013;21(3):450-61. doi:10.1016/j.joca.2012.11.016.

12.

Zhao

Ito

Nakahata

Tai

Saito

, et al. One session of 20 N cyclic compression induces chronic knee osteoarthritis in rats: a long-term study. Osteoarthr Cartil Open. 2022;4:100325. doi:10.1016/j.ocarto.2022.100325.

13.

Ito

Nakahata

Nishitani

Kuroki

Aoyama

. Effects of in vivo cyclic compressive loading on the distribution of local Col2 and superficial lubricin in rat knee cartilage. J Orthop Res. 2021;39(3):543-52. doi:10.1002/jor.24812.

14.

Nakahata

Zhao

Kuroki

Aoyama

Ito

. A non-invasive method for generating the cyclic loading-induced intra-articular cartilage lesion model of the rat knee. J Vis Exp. 2021;(173):10.3791/62660. doi:10.3791/62660.

15.

Sluijs

Geesink

Linden

Bulstra

Kuijer

Drukker

. The reliability of the Mankin score for osteoarthritis. J Orthop Res. 1992;10(1):58-61. doi:10.1002/jor.1100100107.

16.

Love

Huber

Anders

. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550. doi:10.1186/s13059-014-0550-8.

17.

Szklarczyk

Gable

Lyon

Junge

Wyder

Huerta-Cepas

, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47:D607-D613. doi:10.1093/nar/gky1131.

18.

Shannon

Markiel

Ozier

Baliga

Wang

Ramage

, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498-504. doi:10.1101/gr.1239303.

19.

Taminau

Meganck

Lazar

Steenhoff

Coletta

Molter

, et al. Unlocking the potential of publicly available microarray data using inSilicoDb and inSilicoMerging R/Bioconductor packages. BMC Bioinformatics. 2012;13:335. doi:10.1186/1471-2105-13-335.

20.

Johnson

Rabinovic

. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007;8(1):118-27. doi:10.1093/biostatistics/kxj037.

21.

Appleton

CTG

Pitelka

Henry

Beier

. Global analyses of gene expression in early experimental osteoarthritis. Arthritis Rheum. 2007;56(6):1854-68. doi:10.1002/art.22711.

22.

Nam

Perera

Liu

Rath

Deschner

Gassner

, et al. Sequential alterations in catabolic and anabolic gene expression parallel pathological changes during progression of monoiodoacetate-induced arthritis. PLoS ONE. 2011;6(9):e24320. doi:10.1371/journal.pone.0024320.

23.

Sato

Konomi

Yamasaki

Aratani

Tsuchimochi

Yokouchi

, et al. Comparative analysis of gene expression profiles in intact and damaged regions of human osteoarthritic cartilage. Arthritis Rheum. 2006;54(3):808-17. doi:10.1002/art.21638.

24.

Barter

Swan

Rankin

Rowan

Santibanez-Koref

, et al. Identification of the pathogenic pathways in osteoarthritic hip cartilage: commonality and discord between hip and knee OA. Osteoarthritis Cartilage. 2012;20(9):1029-38. doi:10.1016/j.joca.2012.05.006.

25.

Grässel

. The role of peripheral nerve fibers and their neurotransmitters in cartilage and bone physiology and pathophysiology. Arthritis Res Ther. 2014;16(6):485. doi:10.1186/s13075-014-0485-1.

26.

Burgkart

Glaser

Hyhlik-Dürr

Englmeier

Reiser

Eckstein

. Magnetic resonance imaging–based assessment of cartilage loss in severe osteoarthritis: accuracy, precision, and diagnostic value. Arthritis Rheum. 2001;44(9):2072-7. doi:10.1002/1529-0131(200109)44:9<2072::AID-ART357>3.0.CO;2-3.

27.

Salhotra

Shah

Levi

Longaker

. Mechanisms of bone development and repair. Nat Rev Mol Cell Biol. 2020;21:696-711. doi:10.1038/s41580-020-00279-w.

28.

Yik

Acharya

Kumari

Fierro

Haudenschild

, et al. The oncogene LRF stimulates proliferation of mesenchymal stem cells and inhibits their chondrogenic differentiation. Cartilage. 2013;4(4):329-38. doi:10.1177/1947603513497570.

29.

Novack

. Role of NF-κB in the skeleton. Cell Res. 2011;21(1):169-82. doi:10.1038/cr.2010.159.

30.

Kenneth

Miguel

Eleonora

Rosa

Mary

. NF-κB signaling: multiple angles to target OA. Curr Drug Targets. 2010;11:599-613. doi:10.2174/138945010791011938.

31.

Hiraoka

Hosoda

Chiba

Kishimoto

. SGK phosphorylates Cdc25 and Myt1 to trigger cyclin B-Cdk1 activation at the meiotic G2/M transition. J Cell Biol. 2019;218:3597-611. doi:10.1083/jcb.201812122.

32.

Zhong

Yang

Wang

Zheng

, et al. KCTD12 promotes tumorigenesis by facilitating CDC25B/CDK1/Aurora A-dependent G2/M transition. Oncogene. 2017;36:6177-89. doi:10.1038/onc.2017.287.

33.

Huang

Wang

Meng

Liu

Zhao

, et al. Lethal (3) malignant brain tumor-like 2 (L3MBTL2) protein protects against kidney injury by inhibiting the DNA damage-p53-apoptosis pathway in renal tubular cells. Kidney Int. 2018;93:855-70. doi:10.1016/j.kint.2017.09.030.

34.

Liu

Zhou

Sang

Zhang

Zhao

, et al. Silica nanoparticles induce spermatogenesis disorders via L3MBTL2-DNA damage-p53 apoptosis and RNF8-ubH2A/ubH2B pathway in mice. Environ Pollut. 2020;265:114974. doi:10.1016/j.envpol.2020.114974.

35.

Zhang

Lin

Gao

Feng

, et al. Silencing of functional p53 attenuates NAFLD by promoting HMGB1-related autophagy induction. Hepatol Int. 2020;14:828-41. doi:10.1007/s12072-020-10068-4.

36.

Kühn

D’Lima

Hashimoto

Lotz

. Cell death in cartilage. Osteoarthritis Cartilage. 2004;12:1-16. doi:10.1016/j.joca.2003.09.015.

37.

Ito

Maruyama

Sakai

Izumi

Moriishi

Yoshida

, et al. Overexpression of Cdk6 and Ccnd1 in chondrocytes inhibited chondrocyte maturation and caused p53-dependent apoptosis without enhancing proliferation. Oncogene. 2014;33:1862-71. doi:10.1038/onc.2013.130.

38.

Dong

Xia

Wang

Zhao

Tian

. Investigation of candidate genes for osteoarthritis based on gene expression profiles. Acta Orthop Traumatol Turc. 2016;50:686-90. doi:10.1016/j.aott.2016.04.002.

39.

Sun

Yan

Yin

Zhang

. Identification of genes associated with osteoarthritis by microarray analysis. Mol Med Rep. 2015;12:5211-6. doi:10.3892/mmr.2015.4048.

40.

Dirks

Viemann

Beyersdorf

Härtel

Morbach

. Insights into B-cell ontogeny inferred from human immunology. Eur J Immunol. 2023;53(6):e2250116. doi:10.1002/eji.202250116.

41.

Xie

Doody

Shuweihdi

Conaghan

Ponchel

. B-cell capacity for expansion and differentiation into plasma cells are altered in osteoarthritis. Osteoarthritis Cartilage. 2023;31(9):1176-88. doi:10.1016/j.joca.2023.03.017.

42.

Keshvari

Caruso

Teakle

Batoon

Sehgal

Patkar

, et al. CSF1R-dependent macrophages control postnatal somatic growth and organ maturation. PLoS Genet. 2021;17(6):e1009605. doi:10.1371/journal.pgen.1009605.

43.

Wittrant

Gorin

Mohan

Wagner

Abboud-Werner

. Colony-stimulating factor-1 (CSF-1) directly inhibits receptor activator of nuclear factor-κB ligand (RANKL) expression by osteoblasts. Endocrinology. 2009;150:4977-88. doi:10.1210/en.2009-0248.

44.

Garcia

Hartkamp

Malvar-Fernandez

IEV

Lin

Wong

, et al. Colony-stimulating factor (CSF) 1 receptor blockade reduces inflammation in human and murine models of rheumatoid arthritis. Arthritis Res Ther. 2016;18:75. doi:10.1186/s13075-016-0973-6.

45.

Toh

Bonnefoy

Accart

Cochin

Pohle

Haegel

, et al. Bone- and cartilage-protective effects of a monoclonal antibody against colony-stimulating factor 1 receptor in experimental arthritis. Arthritis Rheumatol. 2014;66(11):2989-3000. doi:10.1002/art.38624.

46.

Straub

. Chapter I—history of immunology research. In: Straub

editor. The origin of chronic inflammatory systemic diseases and their sequelae. San Diego (CA): Academic Press; 2015. p. 1-58. doi:10.1016/B978-0-12-803321-0.00001-X.

47.

Malemud

. Negative regulators of JAK/STAT signaling in rheumatoid arthritis and osteoarthritis. Int J Mol Sci. 2017;18:484. doi:10.3390/ijms18030484.

48.

Kleppe

Lahortiga

El Chaar

Keersmaecker

Mentens

Graux

, et al. Deletion of the protein tyrosine phosphatase gene PTPN2 in T-cell acute lymphoblastic leukemia. Nat Genet. 2010;42(6):530-5. doi:10.1038/ng.587.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.09 MB

0.07 MB