Characterization of the single cell landscape in normal and osteoarthritic equine joints

Introduction

Osteoarthritis (OA) represents a major source of disability, affecting 33 million people in the United States and 655 million people worldwide, making it the most prevalent type of arthritis (1,2). In addition to the morbidity associated with OA, the disease has a substantial economic impact, with an estimated cost of $136 billion in healthcare annually (2). There is a high lifetime likelihood of developing symptomatic knee OA in humans (45%), which highlights the need for an improved understanding of OA pathogenesis to facilitate the development of novel diagnostic biomarkers and disease-modifying drugs (3). Similarly, OA is one of the most common disorders treated in equine clinical practice, estimated to affect one third of thoroughbred racehorses and 80% of horses over 15 years of age (4,5). The comparable disease prevalence, similarity in joint volume, cartilage thickness, and articular cartilage loading forces between humans and horses makes spontaneously occurring OA in horses a valuable animal model (6-9). Despite the high prevalence of OA across species, early detection and therapeutic treatment options are limited, due in part to an incomplete understanding of the mechanisms underlying OA progression (10).

The current understanding of OA pathogenesis suggests the innate immune system, particularly myeloid cells, play an important role in regulating and perpetuating low-grade inflammation (11-13). Macrophages are the most common myeloid cell type in synovium and are thought to promote sustained inflammation through release of inflammatory mediators in response to cartilage-derived damage associated molecular patterns (DAMPs) (11-13). The clinical impact of synovial macrophages in OA progression has been demonstrated through ex vivo and in vivo depletion studies (11,14,15). These investigations revealed that a reduction of synovial macrophages resulted in decreased synovial inflammation in experimentally induced OA, emphasizing the key role of macrophages in OA. In addition to macrophages, various T cell subsets have been implicated in OA progression in human, murine, and most recently equine models. These connections between OA pathogenesis and T cell perturbations have been made in the peripheral blood, synovial tissue, and synovial fluid of OA-affected individuals, with increases in abundance and cytokine profile associated with T helper 1 (Th1), Th9, and Th17 cells in disease (16). Recently, evaluation of synovial fluid from equine patients with naturally occurring post-traumatic OA (PTOA) via flow cytometry and enzyme-linked immunosorbent assay revealed an imbalance in the ratio of regulatory T cells and Th17 cells with an increase in Th17 cells in OA (17). Collectively, there is evidence supporting the involvement of multiple immune cell types in OA pathogenesis, but a comprehensive analysis investigating multiple cell types at once is lacking.

The heterogenous inflammatory responses associated with degenerative joint disease pose a significant challenge in the identification of appropriate therapeutic approaches to treat OA (18,19). Single-cell RNA sequencing (scRNA-seq) offers a powerful approach to study the transcriptomic complexity of heterogenous tissues and enables the high-resolution analysis of individual cells. The utility of scRNA-seq has been recently reviewed in the context of orthopedic research (20,21) and has gained significant interest in the past several years in the study of OA (22-25). To highlight a few contemporary studies, Ji et al. recently performed scRNA-seq on chondrocytes from 10 human patients with OA undergoing knee arthroplasty (22), revealing seven transcriptomically distinct chondrocyte populations, including three novel phenotypes predicted to have distinct functions. Mizoguchi et al. [2018] and Zhang et al. [2019] further characterized synovium from rheumatoid arthritis patients using scRNA-seq to characterize fibroblast populations (26,27), and Lewis et al. built on this to classify transcriptional subunits of rheumatoid arthritis synovium into three pathotypes (fibroblastic, macrophage-rich and lymphomyeloid) (28). Tan et al. [2023] implemented scRNA-seq to investigate immune-related biomarkers in synovium of OA patients, identifying vascular endothelial growth factor (VEGF) as a potentially promising biomarker of OA (23). In a mouse model of PTOA with anterior cruciate ligament (ACL) tear, scRNA-seq transcriptomic analyses revealed the most dramatic changes following injury in the monocyte and macrophage cell populations, highlighting alterations in the immune microenvironment (24). These studies in humans and rodent species have generated key data regarding responses in synovium, but there remains a need to investigate how individual cells within synovial fluid are impacted by disease.

Therefore, in this study we applied scRNA-seq to address two objectives. First, we applied the dataset to define the transcriptional profiles of cells within equine synovial fluid and synovium, and second, we investigated how the cellular composition and transcriptomic profiles are altered in the presence of OA. Our analysis revealed the composition of synovial fluid to be entirely immune cells (CD45⁺), with 8 major cell types identified. Within the synovial fluid, we observed an overrepresentation of IL23R⁺ gamma-delta (γδ) T cells enriched in gene expression patterns suggestive of Th17 mediated immunity. Through analysis of macrophages in synovial tissue and fluid, we identified a reduction in the relative abundance of synovial macrophages with a corresponding increase in the relative abundance of synovial fluid macrophages. Overall, this work provides key insights into synovial responses to OA while also indicating myeloid cells and IL23R⁺ γδ T cells are important contributors to disease. We present this article in accordance with the ARRIVE reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-24-40/rc).

Methods

Study overview

Four systemically healthy skeletally mature male castrated Quarter Horses (aged 5, 13, 16, and 16 years), were used for the study (Figure 1, Table 1). Horses of the same breed were enrolled to mitigate potential genetic variability which could introduce biological noise to the dataset. Paired synovial fluid samples (normal tarsus and OA carpus) were obtained from two horses (Horse 1 and 2) and an additional OA sample (carpus) and normal sample (tarsus) were obtained from two separate horses (Horse 3 and 4). Paired synovium samples (normal carpus and OA tarsus) were obtained from one of the horses (Horse 2) to complete a comparison to synovial fluid from the same joints. The tarsus was used as the control joint due to bilateral carpal OA in enrolled cases. Horses were determined to be systemically healthy by physical examination and bloodwork (complete blood count and chemistry panel), and to be free of Salmonella spp. by fecal polymerase chain reaction (PCR). This study was approved by the Institutional Animal Care and Use Committee at Colorado State University (IACUC protocol #926) and conducted according to the national guidelines under which the institution operates and National Institutes of Health (NIH) Guidelines for the Care and Use of Laboratory Animals (8^th edition).

Figure 1 Study overview. scRNA-seq was performed on 3 synovial fluid samples and 1 synovium sample from normal (tarsus) and osteoarthritic joints (carpus) for a total of 8 samples. OA, osteoarthritis; scRNA-seq, single-cell RNA sequencing.

Horse enrollment

Horses were screened for inclusion by physical examination, lameness evaluation performed by two board-certified veterinary surgeons (L.M.P., L.G.), who were aware of which horses and joints were diagnosed with OA and therefore sampled. Radiographs (four-view per limb) of the affected joints (carpus in all cases) were obtained and OA was diagnosed via a combination of clinical examination, lameness/gait evaluation, and radiographic evidence. When possible, a macroscopic scoring system was also used to describe gross changes in the middle carpal joint (29). A traumatic origin was suspected as the triggering factor for all OA cases in this case series. Horses were required to be sound in the unaffected limb and free of radiographic evidence of tarsal OA prior to study enrollment for acquisition of normal joint samples, and to be lame (minimum grade 3/5, American Association of Equine Practitioners’ lameness scoring system) with radiographic evidence of OA in one middle carpal joint for inclusion of the joint as an osteoarthritic sample (30). As horses were enrolled based on inclusion/exclusion criteria for OA diagnosis and not assigned to receive treatment allocation, randomization of experimental units (i.e., horses) was not performed. Additionally, as treatments were not administered, confounders (e.g., order in which treatments were administered) were not controlled.

Sample collection and processing

Synovial fluid (5 to 10 mL) was sampled from the affected middle carpal joint (OA sample) and tibiotarsal joint (normal control sample) of each horse (i.e., experimental unit). Synoviocenteses to obtain synovial fluid were performed implementing aseptic technique using sterile gloves, needles (20 gauge, 1.5 inches), and syringes (31). To obtain a single-cell suspension synovial fluid samples were treated with 30 µg/mL hyaluronidase (Sigma-Aldrich; St. Louis, MO, USA) for 15 min at 37 ℃, washed with phosphate-buffered saline (PBS), then filtered through a 40 µM cell strainer. Red blood cells were lysed using ammonium-chloride potassium (ACK) lysis buffer then cell suspensions were resuspended in PBS with 0.04% molecular grade bovine serum albumin (BSA) (Sigma-Aldrich).

Synovium was obtained in one case (Horse 2) immediately following euthanasia. Briefly, the horse was euthanized with pentobarbital (1 mL/5 kg body weight) administered via indwelling catheter, then the middle carpal and tibiotarsal joints were dissected aseptically and synovium was sampled from four quadrants of each joint (dorsomedial, dorsolateral, palmar- or plantaromedial, palmar- or plantarolateral) and pooled for analysis. Synovium was disaggregated into single cell suspensions as previously described (27). Briefly, synovium was mechanically minced then digested with 200 U/mL Collagenase Type II (Sigma-Aldrich) for 45 minutes at 37 ℃ with agitation every 10 minutes. Digested suspensions were washed with PBS, filtered through a 40 µM cell strainer, then ACK lysis buffer was added for 5 minutes at room temperature to lyse red blood cells. Lastly, a final wash at 100 rcf for 15 minutes was completed to remove small debris and platelets. The prepared cell suspensions were then resuspended in PBS with 0.04% non-acetylated BSA (Sigma-Aldrich).

scRNA-seq processing

The cell concentration for each tissue source was determined using a hemocytometer and trypan blue staining then adjusted to obtain a concentration between 700 to 1,200 cells/µL. Once in solution, cells were processed on a Chromium iX instrument within 2 hours. Briefly, single cells were isolated and tagged with unique cell barcodes using a Chromium iX instrument with a target of 5,000 cells per sample (Next GEM Single Cell 3' Reagent Kit v3.1) (10x Genomics; Pleasanton, CA, USA). The cDNA was then pooled, PCR amplified, cleaned up, and Illumina sequencing adapters were added following manufacture recommended protocols. Sample quality of each cDNA library was evaluated using a LabChip Bioanalyzer (PerkinElmer; Waltham, MA, USA) then submitted for sequencing on an Illumina NovaSeq 6000 sequencer with a target of 100,000 (synovial fluid samples Normal_1 and OA_1) or 50,000 (all other samples) 150 bp paired-end reads per cell (Novogene; Sacramento, CA, USA).

Data pre-processing

Raw sequencing data were demultiplexed, adaptors trimmed, and aligned to the equine genome (EquCab3.0; Ensembl) using a Cell Ranger analysis pipeline with default settings (10x Genomics; version 6.1.2) (Table S1). Cell by gene count matrices were then imported into Seurat (version 4.3.0) to complete dimension reduction and downstream analysis of the generated single-cell gene expression data (32). Briefly, each sample was filtered to retain only cells that fell within quality control parameters for number of transcripts [unique molecular identifier (UMI)] per cell (200< UMI <30,000), number of features per cell (100< nFeature <4,000), and percent of reads that mapped to mitochondrial genes (percent.MT <12.5%). Further filtering was completed using DoubletFinder (version 2.0.3) to identify and remove putative cell doublets (33). Following doublet removal, each sample was log normalized, highly variable genes were identified, and non-linear dimension reduction was completed. All samples (keeping tissue types separate, i.e., synovium and synovial fluid) were then normalized using a SCTransform protocol (Pearson residuals of regularized negative binomial regression) and the data were scaled using “percent.MT” as latent variable in a linear regression framework to minimize the impact of mitochondrial reads on dimension reduction and integration. Canonical correlation analysis (CCA) using conserved variable features between samples was then applied to complete pairwise integration and project the data in the same embedding (34). Unsupervised clustering was completed on each dataset and numerical cluster identifiers were generated such that the cluster containing the greatest number of cells was labeled “0”, and the cluster size decreases incrementally as the cluster identifier increases in value. The size of each cluster within each sample is presented in Table S2 (synovial fluid) and Table S3 (synovial tissue). For all datasets, a Uniform Manifold Approximation and Projection (UMAP) embedding was used to complete dimension reduction. Subset analysis on major cell types was completed through subcluster analysis which consisted of subsetting on the cells of interest, separating samples into biological replicates, then re-integrating the data and repeating unsupervised clustering and dimension reduction. Clustree (version 0.4.4) was used to evaluate cluster stability and determine ideal clustering resolution (35). The key clustering parameters used for each data subset are defined in Table S4.

Cell type annotation

After the data were pre-processed, manual and reference-based methods (singleR; version 1.8.1 and Seurat) were utilized to classify major cell populations (32,36). In addition to human and murine cell type references, an annotated equine peripheral blood mononuclear cell (PBMC) scRNA-seq reference dataset was used to assist in the classification of T cell subtypes (37). The equine PBMC reference dataset was used to complete label transfer using Seurat’s reference mapping workflow (32). Manual annotation of the dataset consisted of evaluation of unsupervised clustering to determine if the cell type subdivisions represented a biologically relevant entity. If substantial overlap in cell type gene signatures were observed, we presumed a portion of the data was over clustered. To address this, clusters were collapsed into populations that more closely fit into biologically relevant cell type divisions. The collapsing of cell populations was informed by clustree. The unsupervised clustering (supplemental figures) and annotated data (primary figures) are presented throughout.

Cell type subsets that did not fit into previously described cell types were annotated using a defining feature (i.e., GZMA, IL23R) followed by well-defined cell subtype name (i.e., γδ T cell, macrophage). Cell type gene lists, as presented in supplemental data, were generated using the FindAllMarkers function (test.use = “wilcox”, logfc.threshold =0.1, min.pct =0.25, only.pos = TRUE, return.thresh =0.01) on each of the major cell type subpopulations (available online: https://cdn.amegroups.cn/static/public/atm-24-40-1.csv) and the 26 final cell types identified throughout analysis (available online: https://cdn.amegroups.cn/static/public/atm-24-40-2.csv). In addition to full cell type gene lists, short gene lists that define each cell type were compiled using a selection of the top features that define each cell type with preference given to unique features in the top 50 genes (weighted by adjusted P value) identified to define each cell population.

Differential abundance analysis

Following cell type identification, shifts in cell type abundances were evaluated using a Monte Carlo permutation test implemented through the R package scProportionTest (version 0.0.0.9000) (38). To complete Monte Carlo permutation testing biological replicates were downsampled to obtain equal representation from each sample, then a null distribution was generated by pooling cells and randomly segregating the cells into clusters 1,000 times. For each permutation, the log₂(fold change) was calculated (contrasting OA to normal within each cluster), then the P value was calculated as the number of permutations that were more extreme than the actual log₂(fold change). P values were then corrected using the false discovery rate (FDR) correction method. An abundance change was considered statistically significant if the actual log₂(fold change) was greater than 0.58 and the FDR corrected P value was less than 0.01. When reporting differential abundance in the main text, the mean percentage ± standard deviation is provided in parentheses.

To further evaluate differential abundances on a by-sample level we used natural log transformed odds ratio (logOR) and Fisher’s exact test to provide evidence that one (or more) sample(s) were over (or under) abundant relative to random distribution across all samples, an approach adopted from Zheng et al. (39). To complete the analysis, we constructed 2×2 contingency tables for each combination of cluster i and sample j. The constructed tables contained the number of cells of cluster i in sample j, the number of cells of cluster i in other samples, the number of cells of non-i clusters in sample j, the number of cells of non-i clusters in other samples. A Fisher’s exact test was then applied to each table and P values were adjusted using the FDR method. ORs were then log transformed to obtain symmetrical distributions and abundance changes were considered statistically significant if the adjusted P value was below 0.01 and |logOR| >0.69 (the equivalent of OR >2 or OR <0.5).

Differential gene expression (DGE) analysis

DGE analysis was investigated by completing pseudobulk conversion followed by a DEseq2 workflow (40). This approach consisted of removing features that had less than 10 cells with a non-zero expression value then collapsing raw counts data for each sample into one column. Each sample was required to be derived from a minimum of 5 cells to be included in the pseudobulk count matrix. If a sample did not meet this threshold, then the sample was excluded from DGE analysis. Statistical testing was conducted differently depending on the application within the study. For analysis comparing gene expression between two cell type clusters the P values were determined by testing the null hypothesis that the |log₂(fold change)| was less than 1. Features were then considered to be significantly differentially expressed if the adjusted (FDR) P value was less than 0.01. For analysis comparing cells from OA tissues to normal tissues each P value was determined by testing the null hypothesis that the |log₂(fold change)| was equal to 0. Features were then considered to be significantly differentially expressed if the adjusted (FDR) P value was less than 0.1 and an |log₂(fold change)| greater than 1.

Any subsequent gene set enrichment analysis (GSEA) was completed using lists of upregulated or downregulated genes and the enricher function from clusterProfiler (version 4.2.2) was used with Reactome or Gene Ontology gene set terms (41,42). Gene sets with an adjusted P value of 0.05 or lower (FDR correction method) were considered significantly enriched. In addition to GSEA, module scoring was completed using the AddModuleScore function from Seurat to evaluate enrichment of gene lists within individual cells.

Cell-cell interaction inference analysis

The R package CellChat (version 1.5.0) was used to make inferences about cell-cell interactions between cell types in synovial tissue and fluid (43). Using a list of known receptor-ligand (R-L) pairs, we calculated the interaction scores (strength and weight) which represent the potential of two cells interacting. Equine gene symbols were converted to human orthologs (only retaining 1:1 orthologs) then analysis was completed using the final annotated datasets to evaluate the interactivity between ligand expressing cells (senders/outgoing signals) and receptor expressing cells (receivers/incoming signals). Inferences regarding potential interactivity were made based on the law of mass action using the average expression values of receptors and ligands. Statistical enrichment of interaction networks within cell types from normal and OA tissues was determined using permutation testing to generate a null distribution based on random sampling. As recommended in the CellChat manual, P values (unadjusted) were calculated then signaling networks with a P value <0.05 were considered significantly enriched.

Statistical analysis

Raw data collected through completion of scRNA-seq on 6 equine scRNA-seq synovial fluid samples and 2 equine synovial tissue samples were generated in this study. Sample collection was completed in batches with normal and OA samples processed side-by-side when possible (exception being OA_3 and Normal_3). The batches were (I) Normal_1 + OA_1; (II) Normal_2 + OA_2; (III) OA_3; and (IV) Normal_3. Batch 2 consisted of both synovial fluid and synovium samples. The batch correction protocol for data integration, dimension reduction, and clustering treated each sample as its own batch. A blocking design was used when completing DGE analysis with DEseq2 in which we paired OA_1 and Normal_1, OA_2 and Normal_2, and treated OA_3 and Normal_3 as separate batches. Biological replicates were used for pseudobulk DGE analysis, while cellular replicates were used for all other analysis completed in this study. Detailed descriptions of the statistical analyses and significance thresholds used in this study are provided in the respective methods section.

Results

scRNA-seq of equine synovial fluid reveals an overrepresentation of myeloid cells in horses with OA

Analysis of 28,639 cells from the synovial fluid of three normal and three osteoarthritic joints revealed the presence of 8 major cell types which included: CD8 T cells, γδ T cells, CD4 T cells, macrophages, dendritic cells (DCs), neutrophils, B cells and cycling cells (Figure 2A,2B, available online: https://cdn.amegroups.cn/static/public/atm-24-40-1.csv). Briefly, T cells were defined by the expression of CD3E plus GZMA and CTSW (CD8 T cells), TRDC and KLRB1 (γδ T cells), or diffuse expression of CD4 coupled with a lack of CD8 and γδ T cell associated features (CD4 T cells). Macrophages and DCs were identified as AIF1 (Iba1) expressing cells with macrophages additionally expressing CD68 and DCs defined by FLT3/CD1C expression. Neutrophils were classified based on the lack of DRA (MHC class II) expression with high expression of CSF3R and S100A12. B cells were annotated based on MS4A1 (CD20) expression. Lastly, a heterogenous population of cells with gene expression patterns suggestive of actively going through the cell cycle (TOP2A, CENPF, PCNA) was identified (Cycling cells).

Figure 2 Single-cell RNA sequencing reveals 8 major cell type populations in equine synovial fluid. (A) UMAP representation of the 28,639 cells obtained from equine synovial fluid colorized by assigned major cell type identity. (B) Expression of canonical cell type features used in cell classification. (C) UMAP embedding colorized by sample from which the cell came from. Purple colors indicate cells came from osteoarthritic synovial fluid while green colors indicate cells came from normal synovial fluid. (D) Heatmap depicting the number of DEGs comparing cells from OA to normal synovial fluid within each cell type. Features are included in counts if adjusted P value <0.1 and |log₂(fold change)| >1. Features with a positive log₂(fold change) are in right column (“Up”), while features with a negative log₂(fold change) are in left column (“Down”). (E) Volcano plot depicting results of pseudobulk differential gene expression analysis comparing OA to normal cells within the cluster of actively cycling cells. Top 20 genes (weighted by adjusted P value) are labeled in red (“Up in OA”) and blue (“Up in Normal”). (F) Dot plot depicting the results of module scoring for non-cycling cell types (y-axis). The x-axis terms, “Cycling_down” and “Cycling_up” correspond to the DEGs identified in (E) labeled as “Up in Normal” and “Up in OA”, respectively. Throughout the figure “gd” T cell refers to γδ T cell. UMAP, Uniform Manifold Approximation and Projection; OA, osteoarthritis; DEG, differentially expressed gene; γδ, gamma-delta.

Following annotation of major cell types, we evaluated OA-associated changes in cell type abundances. We observed increased relative abundances of macrophage (normal =16.7%±8.4%; OA =37.2%±16.0%) and DCs (normal =9.36%±5.4%; OA =19.2%±7.8%) as well as reductions in the abundances of CD4 T (normal =20.3%±6.7%; OA =8.3%±5.3%) and CD8 T (normal =32.4%±6.5%; OA =14.5%±8.1%) in OA synovial fluid relative to normal (Figure 2C, Figure S1, Table S5). We next completed DGE analysis between cells from OA and normal joints within each of the 8 major populations (available online: https://cdn.amegroups.cn/static/public/atm-24-40-3.csv). This analysis demonstrated that the population of cycling cells exhibited the greatest number of differentially expressed genes (DEGs), whereas lymphocytes (B/T cells) were found to have minimal detectable DEGs (Figure 2D). Within the cycling cells we identified myeloid associated features (AIF1, TIMP1, DRA, CST3) to be enriched and T cell associated features (CTSW, CCL5, and GZMA) to be downregulated in OA joints (Figure 2E). Because gene expression patterns and cell type abundance shifts are closely related, we next used the genes found to be upregulated in OA (“Cycling_up”) and upregulated in normal (“Cycling_down”) to complete module scoring of the 7 non-cycling major cell types. This analysis indicated that the genes enriched in OA within the cycling cell population mapped to neutrophils, DCs, and macrophage, whereas the downregulated (enriched in normal cells) genes mapped to CD8, CD4, and γδ T cells (Figure 2F). This enrichment pattern suggests myeloid cells make up a larger relative proportion of cycling cells in osteoarthritic synovial fluid. In summary, we defined the cellular heterogeneity within equine synovial fluid and identified a population shift in the OA joint which favors a myeloid cell predominance. Next, we further investigated heterogeneity within the major cell types through completion of subcluster analysis.

Myeloid cell transcriptional programs are altered in OA

Higher resolution analysis of myeloid cells within equine synovial fluid using subcluster analysis provided evidence of 12 transcriptionally distinct cell types (Figure 3A,3B, available online: https://cdn.amegroups.cn/static/public/atm-24-40-4.csv). We identified 4 macrophage populations, which were classified as GPNMB⁺ macrophages, activated (CD5L⁺) macrophages, CCL2⁺ macrophages, and TPPP3⁺ macrophages. In addition to macrophages, a population of monocytes (CXCL8⁺/IL18⁺) and intermediate cells with transcriptional programs suggestive of transitioning from monocytes to macrophages (Mo-Mac) were identified. Evaluation of the myeloid populations in the context of traditional dichotomous pro-/anti-inflammatory nomenclature indicated that monocytes exhibited gene expression patterns consistent with a pro-inflammatory profile, while GPNMB⁺ macrophage had gene expression patterns most consistent with anti-inflammatory macrophage (Figure S2A) (44). The 5 DC subsets identified in the dataset were annotated as conventional type 1 DCs (cDC1; FLT3⁺/DNASE1L3⁺/BATF3⁺), conventional type 2 DCs (cDC2; FLT3⁺/CD1C⁺), plasmacytoid DCs (pDC; FLT3⁺/TCF4⁺/MS4A1⁺), migratory DCs (migDCs; FLT3⁺/LAMP3⁺/CCR7⁺) and cycling DCs (FLT3⁺/TOP2A⁺). Additionally, a cluster of neutrophils (DRA^low, CSF3R⁺, SELL⁺) was identified.

Figure 3 Unsupervised clustering reveals 12 distinct macrophage clusters in equine synovial fluid. (A) UMAP representation of 13,173 myeloid cells colorized by unique cluster identity as determined through unsupervised clustering. (B) Dot plots depicting the expression of features enriched in the different macrophage subpopulations. (C) Differential abundance analysis of macrophage subpopulations in osteoarthritic relative to normal synovial fluid (y-axis indicates the log₂(FC) in cell type abundance and x-axis is the cell type being evaluated. Significance was determined using Monte Carlo permutation testing with adjusted P value threshold set to 0.01 and |log₂(FC)| >0.58. The lower heatmap depicts the natural logOR with significantly under or overrepresented cell types (relative to a random distribution across all samples) indicated with an *. Statistical significance was classified based on an adjusted P value <0.01 and |logOR| >0.69 (the equivalent of OR >2 or OR <0.5). (D) Volcano plot depicting differential gene expression between OA and normal synovial fluid within all myeloid cells. The significantly DEGs are labeled in red (“Up in OA”) and blue (“Up in Normal”). (E) Heatmap depicting the number of DEGs between OA and normal synovial fluid within each cell type. Features are included in counts if adjusted P value <0.1 and | log₂(FC)| >1. Features with a positive log₂(FC) are in right column (“Up”), while features with a negative log₂(FC) are in left column (“Down”). (F) Heatmap depicting log normalized average expression within each myeloid cluster separated by condition (“OA” or “Normal”). Asterisks indicates statistically significant change in expression between condition adjusted P value <0.1 and |log₂(FC)| >1. The asterisk is located on the condition that had significantly higher expression. (G) UMAP split by cell source (normal = top row, OA = bottom row) that depicts the expression of select upregulated (GPNMB and CPM) and down-regulated (MARCO, SAMD9L) features within all myeloid subpopulations. UMAP, Uniform Manifold Approximation and Projection; Mac; macrophage; FC, fold change; n.s., not significant; logOR, log transformed odds ratio; OA, osteoarthritis; DEG, differentially expressed gene.

Differential abundance analysis using a Monte Carlo permutation test revealed a shift in DC populations in which cDC2s were overrepresented (normal =14.4%±8.7%; OA =23.8%±2.9%) and cDC1s were underrepresented (normal =15.6%±6.1%; OA =6.6%±8.1%) in OA joints relative to normal joints (Figure 3C, Figure S2B, Table S6). When the abundance changes were evaluated at the sample level, we identify contradictory responses within cDC1s from the OA synovial fluid suggesting variability in cDC1 abundances across the 3 OA samples. Specifically, two of the OA samples (OA_2 and OA_3) were driving the reduction, while the third sample (OA_1) was overabundant relative to a random distribution. To further describe equine cDC1s and cDC2s, we completed DGE analysis which provided a complete gene list to facilitate differentiation of the two most abundant DC populations (Figure S2C, available online: https://cdn.amegroups.cn/static/public/atm-24-40-5.csv).

Next, we completed DGE analysis between myeloid cells from OA and normal synovial fluid cells. The analysis detected 6 upregulated and 7 downregulated features (Figure 3D, available online: https://cdn.amegroups.cn/static/public/atm-24-40-6.csv). Subsequent visualization of DEGs in the UMAP embedding revealed features enriched in OA (CD1A3, IL1RN, and ACTA2) to be localized to the region of cDC2s, further supporting overrepresentation of cDC2s in OA (Figure S2D). To minimize the impact of abundance changes on the analysis, we repeated DGE analysis within 9 of the 12 myeloid cell types (cell numbers prohibited analysis of clusters 9–11). Relatively few DEGs were identified, but the analysis identified CD5L⁺ macrophages, CCL2⁺ macrophages, and cDC2s to exhibit the greatest number of DEGs (Figure 3E,3F). In particular the expression of GPNMB and CPM was enriched in OA CD5L⁺ macrophages, while MARCO expression was downregulated in OA cells within cDC2s, CD5L⁺ macrophages, and TPPP3⁺ macrophages (Figure 3G). Together our analysis of synovial fluid myeloid cells suggests that OA impacts the ratio of cDC2 to cDC1s and modulates MARCO, GPNMB, and CPM expression within a subset of myeloid cells.

IL23R expressing γδ T cells are enriched in the synovial fluid of horses with OA

Subcluster analysis of CD4, CD8, and γδ T cells revealed the presence of 15 transcriptomically distinct clusters that were annotated as 12 discrete cell types (Figure 4A, Figure S3A,S3B, available online: https://cdn.amegroups.cn/static/public/atm-24-40-7.csv). Consistent with a previous report using scRNA-seq to describe circulating equine T cells, the division of T cell clusters into CD4 and CD8 was not possible based exclusively on transcript abundance of CD4/CD8A (Figure S3C) (37). To overcome this limitation, we completed reference mapping to an equine PBMC dataset which provided support for the division of T cells into CD4 and CD8 populations (Figure S3D) (37). CD4 T cells present in synovial fluid were most resemblant of “T CD4⁺ non-naïve” cells identified in peripheral blood, while the majority of CD8 T cells were resemblant of “T CD8⁺ memory” (Figure S3D).

Figure 4 IL23R⁺ γδ T cells are overrepresented in osteoarthritis-affected synovial fluid. (A) UMAP representation of 12,653 T cells colorized by assigned cell type identity. (B) Feature plots and (C) dot plots depicting expression of cluster defining features. (D) Differential abundance analysis of T cell subpopulations in osteoarthritic relative to normal synovial fluid (y-axis indicates the log₂(FC) in cell type abundance and x-axis is the cell type being evaluated. Significance was determined using Monte Carlo permutation testing with adjusted P value threshold set to 0.01 and | log₂(FC) | >0.58. The lower heatmap depicts the natural logOR with significantly under or overrepresented cell types (relative to a random distribution across all samples) indicated with an *. Statistical significance was determined based on an adjusted P value <0.01 and |logOR| >0.69 (the equivalent of OR >2 or OR <0.5). Volcano plots depicting DGE analysis between (E) IL23R⁺ γδ T cell 1 and (F) IL23R⁺ γδ T cell 2 versus all other T cells. (G) Gene set enrichment analysis using enriched features from DGE analysis comparing both IL23R⁺ γδ T cell populations to all other T cells. Count (size of dot) indicates the number of features mapping to the gene set. Throughout the figure “gd” T cell refers to γδ T cell. UMAP, Uniform Manifold Approximation and Projection; FC, fold change; n.s., not significant; logOR, log transformed odds ratio; OA, osteoarthritis; γδ, gamma-delta; DGE, differential gene expression.

With CD4 T cells identified, we were able to further divide this population into naïve (SELL⁺/CCR7⁺; c2), activated (CD40L⁺; c1), and regulatory T cells (FOXP3⁺/CTLA4⁺; c6; T_reg) (Figure 4A-4C, Figure S3E). Unexpectedly, the cluster annotated as T_reg demonstrated partial expression of CTLA4/FOXP3, so there may be heterogeneity within the cluster (Figure 4B,4C). Four distinct CD8 T cell populations were identified which included effector T cells (CCL5/GZMA; c0), T cells enriched in interferon gene signatures (ISG15/RSAD2; c11; T-IFN), CX3CR1⁺ CD8 T cells (c5), and DAPL1⁺ CD8 T cells (c3). When completing reference mapping to an equine PBMC reference, CX3CR1⁺ CD8 T cells were found to closely resemble the cell type defined as “PRF1⁺ non-annotated” in Patel et al. (Figure S3D) (37). DAPL1⁺ CD8 T cells represent a poorly demarcated cell type that did not fit into any classical T cell subset and were annotated based on expression of DAPL1, a gene that has been associated with CD8 T cell exhaustion (45).

The last major T cell subset was a population of CD3⁺ cells that clustered distantly from CD4/CD8 T cells and were most consistent with γδ T cells based on T Cell Receptor Delta Constant (TRDC) expression (46). Despite statistical enrichment of TRDC expression, two of the five clusters (c4 and c7) exhibited low TRDC expression relative to the other γδ T cell populations (Figure 4B,4C). As such, c4 was annotated based on the overexpression of GNLY (GNLY_T), while c7 was annotated based on the over expression of KLRD1 (KLRD1_T). Two of the remaining three clusters, c9 and c10, were both defined by the expression of IL23R and SCART1 despite clustering in distinct regions of the UMAP. These two clusters were annotated as IL23R⁺ γδ T cell population 1 (IL23R_gd_T1) and IL23R⁺ γδ T cell population 2 (IL23R_gd_T2), respectively. We completed DGE analysis to investigate the distinguishing features between the two IL23R⁺ γδ T cell populations and identified two features (ENSECAG00000038560 and BLK) to be enriched in IL23R⁺ γδ T cell population 2, but no features were identified to be enriched in IL23R⁺ γδ T cell population 1 (Figure S4A, available online: https://cdn.amegroups.cn/static/public/atm-24-40-5.csv). This finding supports that the two IL23R⁺ γδ T cell populations clustered apart due to the expression of distinct T cell receptor gamma constant regions, as the orthologue of ENSECAG00000038560 is T cell receptor gamma constant 2. The final cluster, c8, was also defined by TRDC, BLK, and ENSECAG00000038560, but lacked IL23R and SCART1 expression. This expression pattern, coupled with overexpression of GZMA, led to the annotation of the cluster as GZMA⁺ γδ T cells (GZMA_gd_T).

Following cell type identification, we completed differential abundance analysis between OA and normal synovial fluid which revealed a marked overrepresentation of IL23R⁺ γδ T cells (IL23R⁺ γδ T1: normal =1.4%±1.2%, OA =4.3%±1.7%; IL23R⁺ γδ T2: normal =0.7%±0.6%, OA =3.2%±2.2%) as well as a more subtle reduction in DAPL1⁺ CD8 T cells (normal =15.6%±6.1%; OA =6.6%±8.1%) and SELL⁺ CD4 T cells (normal =15.6%±6.1%; OA =6.6%±8.1%) (Figure 4D, Figure S4B,S4C, Table S7). Given the detected overrepresentation of IL23R⁺ γδ T in OA synovial fluid, we next sought to better define their transcriptomic signatures, as they may play an important role in OA pathobiology. As such we completed DGE analysis contrasting each IL23R⁺ γδ T cell population with all other T cells (Figure 4E,4F, available online: https://cdn.amegroups.cn/static/public/atm-24-40-5.csv). The analysis revealed 18 overlapping features in the cell type gene signatures, further suggesting the two populations are highly similar (Figure S4D). Features unique to IL23R⁺ γδ T1s included CD70, CD160, and CD27, while features unique to IL23R⁺ γδ T2s were IL2RA, CCR6, TRAT1, and HOPX. To make inferences regarding the inferred cellular function of these two populations we used the gene ontology database to complete GSEA. We found the gene signature for both IL23R⁺ γδ T cell populations to be enriched in cytokine activity, adaptive immune responses, and Th17 differentiation (Figure 4G, Figure S4E,S4F, available online: https://cdn.amegroups.cn/static/public/atm-24-40-8.csv). Th17 cells are thought to play a role in autoimmunity and inflammatory processes (47), so the association with Th17 differentiation coupled with the overrepresentation in OA makes IL23R⁺ γδ T cells an intriguing cell type that warrants further study. In summary, we described the T cell populations present in equine synovial fluid and identified IL23R⁺ γδ T cells (which were defined by gene expression patterns associated with Th17 immunity) to be overrepresented in osteoarthritic joints.

Synovial fluid cells from osteoarthritic joints exhibit increased myeloid cell associated interactions

After identifying 26 transcriptionally distinct cell populations in equine synovial fluid (Figure S5A, Table 2), we next investigated how the presence of OA impacted cell-cell interactions within joint fluid. We used CellChat to make inferences about cell-cell communication activity, and this analysis led to the identification of 1,483 and 1,561 inferred interactions in normal and OA synovial fluid, respectively (Figure S5B) (43). Further investigation of which cell types were driving the increased interactivity in OA revealed that GPNMB⁺ macrophages, monocyte-macrophage intermediate cells (Mo-Mac), and TPPP3⁺ macrophages were the main drivers of the predicted increase in cell-cell interactions (Figure 5A). Evaluation of the interaction strength between cells in joints with OA compared to normal joints suggested the increased interactivity in the aforementioned macrophage populations was largely mediated by myeloid-myeloid cell interactions (mean differential interaction strength =+0.06) and myeloid-lymphocyte signals (mean differential interaction strength =+0.04) (Figure 5B, Figure S5C). Additionally, we observed an overall reduction in lymphocyte-lymphocyte interactions (mean differential interaction strength =−0.02) within OA-affected synovial fluid, which suggests that myeloid interaction networks are more likely to be impacted by OA.

Figure 5 Cell-cell interaction inference analysis reveals increased myeloid associated cell interactions in osteoarthritis-affected synovial fluid. (A) Scatter plot depicting the outgoing (x-axis) and incoming (y-axis) interaction strength for each cell type in normal (top) and osteoarthritic (bottom) synovial fluid. The size of the dot indicates the number of predicted interactions, with larger dots corresponding to more interactions. (B) Heatmap depicting the change in inferred interaction strength between cell types comparing cells in OA to normal synovial fluid. (C) Stacked bar graph illustrating the proportion of interactions contributed by cells from normal (green) and OA (purple) joints for each enriched interaction network. Circos plot depicting the inferred interactions between cell types within enriched networks identified in normal (D; HGF/IGF) and OA (E; FASLG/VISFATIN) synovial fluid. The arrow origin (flat base) represents expression of a ligand, while the arrow terminus (pointed tip) indicates expression of a receptor. Cell types not involved in the network are greyed out. Throughout the figure “gd” T cell refers to γδ T cell. OA, osteoarthritis; HGF, hepatocyte growth factor; IGF, insulin-like growth factor; FASLG, FAS ligand; VISFATIN, nicotinamide phosphoribosyltransferase (aka NAMPT); γδ, gamma-delta.

In total, our analysis predicted 28 receptor-ligand pathways to be active in normal and OA synovial fluid, with 23 of the pathways conserved across the two conditions (Figure 5C, Figure S5D). The top two interaction networks unique to cells in normal synovial fluid were hepatocyte growth factor (HGF) and insulin-like growth factor (IGF) (Figure 5D, Figure S5E). HGF signaling has been associated with increased expression of ARG-1 mRNA and secretion of IL-10 and TGF-β1 (48), while IGF signaling has been associated with programming macrophages to favor oxidative phosphorylation and anti-inflammatory properties (49). Thus, enrichment of these interaction networks within normal synovial fluid suggests the presence of an anti-inflammatory cellular environment in normal joints. The top 3 interaction networks enriched in OA synovial fluid were JAM, VISFATIN, and FASLG (Figure S5F,S5G). Increased FASLG network activity was driven by the upregulation of FAS on GZMA⁺ γδ T cells and IFN-T cells, indicating other T cell populations may be targeting the FAS⁺ cells for induction of apoptosis (Figure 5E, Figure S5F) (50). Enrichment of the VISFATIN network in OA synovial fluid cells was largely mediated by increased expression of ITGA5 and NAMPT (also known as VISFATIN) across myeloid cells and was predicted to signal multiple macrophage, DC, and cycling cell populations (Figure S5G). Previous reports have identified VISFATIN as a pro-inflammatory mediator that plays a role in monocyte differentiation, providing evidence that VISFATIN cell-cell communications may be acting to shape the OA microenvironment (51,52). In summary, our analysis of cell-cell interactions within synovial fluid revealed that myeloid cells are a key cell type involved in the OA-associated increase in cell-cell communication.

scRNA-seq of synovial tissues from one joint with OA compared to tissue from one normal joint reveals a reduction of myeloid cell abundances

To supplement our analysis of synovial fluid we also completed scRNA-seq on synovial tissue from one OA-affected joint (n=10,079 cells) and the synovial tissue from one normal joint (n=7,611 cells). This dataset was largely composed of CD45⁻ synoviocytes, while myeloid cells (primarily macrophages and DCs), T cells, myofibroblasts, endothelial cells, and fibroblasts were also identified (Figure 6A,6B, Table S8, available online: https://cdn.amegroups.cn/static/public/atm-24-40-9.csv). Relative to the synovial fluid, we determined synoviocytes (DNASE1L3, PCOLCE, NECTIN3), fibroblasts (FBLN1, COL1A1, SFRP4), myofibroblasts (ACTA2, NOTCH, MCAM), and endothelial cells (VWF, ESAM1) were only found in synovial tissue (27). The synoviocytes were further subdivided into 6 subclusters which were defined by differential expression of CXCL14, IL13RA2, DNASE1L3, and CD55 (Figure 6A, Figure S6).

Figure 6 Relative abundances of myeloid cells are reduced in osteoarthritis-affected synovium. (A) UMAP representation of 17,541 synovial tissue cells colorized by assigned cell type identity. (B) Feature plots depicting expression of genes enriched in the different synovial tissue cell subsets. (C) Differential abundance analysis comparing osteoarthritic to normal within synovial tissue cells [y-axis indicates the log₂(FC) in cell type abundance and x-axis is the cell type being evaluated]. Significance was determined using Monte Carlo permutation testing with adjusted P value threshold set to 0.01 and |log₂(FC)| >0.58. (D) Heatmap depicting the change in inferred interaction strength between cell types comparing cells in OA to normal synovial tissue. (E) Scatter plot depicting the outgoing (x-axis) and incoming (y-axis) interaction strength for each cell type in normal (left) and OA (right) synovium. The size of the dot indicates the number of predicted interactions, with larger dots corresponding to more interactions. (F) Circos plot depicting the inferred interactions between cell types within enriched networks identified in normal (PDGF/NOTCH) synovial tissue. The arrow origin (flat base) represents expression of a ligand, while arrow terminus (pointed tip) indicates expression of a receptor. Cell types not involved in the network are greyed out. UMAP, Uniform Manifold Approximation and Projection; FC, fold change; n.s., not significant; OA, osteoarthritis; PDGF, platelet-derived growth factor; NOTCH, Notch.

With synovial tissue cell types defined we next evaluated the dataset for OA-associated population shifts. Contrary to what was observed in synovial fluid, the analysis revealed myeloid cell numbers were generally reduced in OA affected synovial tissue (520 cells; 5.2%) relative to normal synovium (1,050 cells; 13.9%) (Figure 6C, Table S9). To further investigate myeloid cell heterogeneity and the cell type similarities between synovial fluid and tissue we integrated myeloid cells from both tissue sources (Figure S7A). The analysis indicated most synovial tissue myeloid cells closely resembled synovial fluid macrophage subpopulations (Figure S7B). Of the macrophage populations present in the synovium, the most abundant population in both normal and OA-affected joints was activated (CD5L⁺) macrophages (42% of synovial myeloid cells) (Figure S7C).

Having identified an increase in predicted cell-cell interactions within synovial fluid, we next applied the same approach to make inferences about alterations in synovial tissue. We identified a similar number of inferred interactions in osteoarthritic (1,094 interactions) and normal (1,075 interactions) synovium, with a reduction in the predicted interaction strength of cells present in OA-affected tissue (Figure 6D,6E, Figure S8A,S8B). The reduced strength of inferred interactions was driven by a reduction in both outgoing (ligand expression) and incoming (receptor expression) within all cell types identified in synovial tissue, suggesting cell-cell cross talk was generally less active in OA joints relative to normal joints (Figure 6E). Despite the broad reduction in inferred interaction strength, synoviocytes and fibroblasts were consistently the cell types with the greatest interaction potential in health and disease. Of the receptor-ligand networks enriched in normal synovium, most were associated with endothelial cell function and vascularization (NOTCH, PDGF, ICAM) (Figure 6F, Figure S8C,S8D). Reduced vascularization in subchondral bone and associated ischemic changes have been associated with OA progression (53,54), which our data provide indirect molecular evidence of vascular changes in the synovium of OA-affected joints.

Discussion

Spontaneously occurring OA in horses is a useful preclinical animal model to investigate the dynamics of synovial immune responses (55). As such, we performed scRNA-seq of synovial fluid and synovium from normal and osteoarthritic joints of horses with naturally occurring disease. Consistent with previous reports using flow cytometry (17), our analysis revealed a heterogeneous synovial fluid cellular microenvironment that was composed entirely of immune cells. Within the synovial fluid cell populations, we identified OA-associated changes which included a relative increase in myeloid cell abundances and an expansion of IL23R expressing γδ T cells. Through integrated analysis with synovium and synovial fluid, we observed that the relative abundance of myeloid cells was reduced in synovium despite a concurrent increase myeloid cell abundances in the synovial fluid of OA joints. Overall, the data presented here provide molecular insights into OA-associated changes within synovium and synovial fluid which provides a valuable resource for future investigations using the equine joint disease model.

Annotation of cell types present in scRNA-seq datasets is often a rate-limiting step in data analysis. This is particularly true for analysis of data obtained from non-traditional animal models due to a lack of species-specific cell type gene signatures. Annotation of our dataset relied on human and murine gene signatures as well as references to equine specific scRNA-seq datasets (27,36,46,56). When possible, annotated equine datasets were used to complete reference mapping or cell type gene signature enrichment scoring to establish connections between previous reports (37). Our final cell type annotations represent identities that are suspected to correspond to biologically relevant subdivisions of cellular populations. Further functional and transcriptomic investigation is required to confirm the cell type identities. To facilitate future application of the annotations presented here, we provide short cell type gene signatures, full gene signatures, and access to the fully annotated dataset for application in label transfer to new equine scRNA-seq datasets (see Data Sharing Statement).

Traditionally, macrophages have been dichotomized into pro-inflammatory and anti-inflammatory groups, but with recent advances there is profound evidence that macrophage phenotypes are more complex than a dichotomy (57,58). The application of a dichotomous annotation approach failed to capture the heterogeneity in our dataset, although it did provide evidence that monocytes represented a pro-inflammatory population and GPNMB⁺ macrophages most likely represented anti-inflammatory macrophages. Of note, the identification of pro-inflammatory monocytes in our analysis, and previous reports of multiple monocyte populations in murine joints, highlights the ability of single-cell transcriptomics to distinguish between infiltrating monocytes and macrophages (24). Outside of macrophage heterogeneity our analysis indicated the presence of an OA-associated shift in DC populations, specifically an increase in cDC2s and a reduction in cDC1s in osteoarthritic joints. Interestingly, a previous report in human OA identified an increase in CD11c⁺ cells in OA-affected joints relative to healthy joints (59). Although the authors discussed this as an enrichment of pro-inflammatory macrophages, our analysis revealed that in horses, both DCs and most macrophage populations express CD11c (ITGAX). As such, our high-resolution analysis provides additional data to suggest that increases in CD11c expressing cells observed in OA-affected joints may indicate more than just a shift in macrophage polarization. To this end, cDC1s are generally described in promoting Th1 immunity while cDC2s are typically described as exhibiting broader functions that can drive Th1, Th2, Th17, and T_reg mediated immunity (60), which may be related to the expansion of Th17-like IL23R⁺ γδ T cells reported in this study.

The balance of T cell mediated immunity in degenerative arthritis and rheumatoid arthritis is key to regulating disease progression. Of specific interest, IL-17 producing T helper (Th17) cells have recently been identified as a cell type that can indirectly promote synoviocyte proliferation and cartilage destruction via inhibition of autophagy (61,62). An earlier study in a rodent model of PTOA demonstrated that senescent OA chondrocytes promote naïve T cell differentiation to Th17 cells which in turn acts to promote a Th17-type immune response (63). Despite the inability to identify a distinct Th17 population in our dataset, we were able to characterize two IL23R⁺ γδ T cell populations. IL23R⁺ γδ T cells have been reported to represent a weakly activated γδ T cell subset which functions as an IL-23 sink to suppress IL-23 mediated Th17 differentiation (64). Therefore, the expansion of two IL23R⁺ γδ T cell populations in equine osteoarthritic joints may represent a compensatory mechanism to quell Th17-mediated inflammatory processes. In support of this phenomenon, neutralizing IL-17 antibodies given in a murine model of collagen-induced arthritis resulted in a reduction of degenerative changes (65). In the context of equine PTOA, Keller et al. [2023] demonstrated an increase in Th17 cells in synovial fluid of horses with severe disease (17). Collectively, there is a growing body of literature supporting the role of Th17 mediated immunity in the OA disease process. Therefore, further exploration of IL-17 or IL-23 targeted therapeutics is warranted in the context of equine OA and with translational relevance to other species suffering from similar disease processes.

The data presented in this study provided the ability to examine cellular and transcriptomic changes in both synovium and synovial fluid. As such, our analysis detected a broad expansion of myeloid cells in synovial fluid, with a corresponding reduction of myeloid cells in the synovium of OA joints. This conflicting finding may represent a breakdown of the immunological barrier present in the joint capsule resulting in migration of myeloid cells from the synovium to synovial fluid, a process that is limited in healthy joints (66,67). While our analysis provides indirect evidence that synovial macrophages migrate into the synovial fluid, it is also possible that infiltrating monocytes could have contributed to this relative increase of myeloid cells in the synovial fluid of OA joints. Further investigation into the potential mechanism of myeloid cell migration is necessary to draw conclusions on the migration hypothesis presented here. Through completion of cell-cell interaction we were able to provide further evidence that myeloid cells play a key role in OA progression. Specifically, we observed a relative increase in the interactivity between myeloid cell populations and other immune cells within the synovial fluid of OA joints. Although the reported increase in interactivity may have been influenced by the cell type abundance shifts associated with OA, the observed increase in predicted interactivity provides evidence that myeloid cells influence cell-cell interactions within the synovial fluid microenvironment. Lastly, in support of the validity of our cell-cell interaction analysis, we identified an increase in vascular signaling network activity within synovium which is consistent with previous reports implicating vascular changes in association with OA (53,54).

The greatest limitation to the current study was the small number of biological replicates included in our dataset, especially in the context of synovium. Furthermore, only male horses were sampled in the study, so any sex differences were not considered in this report. Despite the limited biological replicates, the number of cells obtained in both tissues enabled the high-resolution characterization of cell types found in each tissue. The use of paired samples from the OA-affected and normal limbs from the same horses (in all but one pair of OA-affected and normal samples) contributed to a reduction of batch effects that could have potentially confounded the analysis presented here. Despite this, the methodology used to evaluate changes in relative abundances were exposed to many potential sources of bias (small sample size, mix paired-unpaired study design, and low abundance of rare cell types) that could have impacted the reliability of our findings. Further investigation using larger sample sizes and alternative means of identifying the populations found to be differentially abundant in OA synovial fluid and synovium are necessary to substantiate the findings reported in this study. Furthermore, recent reports have demonstrated that sampling bias can be introduced by completing analysis on distinct joints (68), so the fact that OA samples were collected from the middle carpal joint and the normal samples from the tibiotarsal joint may have impacted the findings. This approach was taken due to the presence bilateral carpal OA in the enrolled cases which prohibited collection from a normal contralateral carpal joint. An additional limitation is that CellChat (cell-cell interaction) analysis relied on a human specific receptor ligand database which may not translate directly to equine cell types. To overcome this limitation, we converted equine gene symbols to human orthologs, but validation of key pathways is an important next step.

Conclusions

The dataset presented here provides annotated scRNA-seq reference datasets of equine synovial fluid and synovium which can be applied in future studies to enhance the utility of the equine OA model. Our comparison of osteoarthritic and normal tissues revealed an expansion of IL23R⁺ γδ T cells and myeloid cells in the synovial fluid of OA-affected joints, which suggests these populations are active in disease progression. Through integrated analysis of synovial fluid and synovium samples, we were able to identify inverse responses in terms of myeloid cell abundances, suggesting that myeloid cells from synovial tissue may migrate into the synovial fluid under diseased conditions. This work provides preliminary evidence that therapeutic strategies targeting IL-17 and IL-23 pathways are worth further investigation in equine OA, with translational relevance to humans and other species.

Acknowledgments

The authors gratefully acknowledge the assistance of Jen Daniels, Ryan Shelton, Natalie Lombard, and other staff members of the Orthopaedic Research Center for their assistance in data collection and care of university-owned research horses.

Funding: This study was funded by the Foundation for the Horse, Grayson Jockey Club Research Foundation, Colorado State University College Research Council Interdisciplinary Pilot Award, NIH/NCATS CTSA 5TL1TR002533-02 (L.M.P.), NIH 5T32OD010437-19 (L.M.P.), NIH/NCATS Colorado CTSA T32TR004366 (Z.W.), and Carolyn Quan and Porter Bennett.

Footnote

Reporting Checklist: The authors have completed the ARRIVE reporting checklist. Available at https://atm.amegroups.com/article/view/10.21037/atm-24-40/rc

Data Sharing Statement: Available at https://atm.amegroups.com/article/view/10.21037/atm-24-40/dss

Peer Review File: Available at https://atm.amegroups.com/article/view/10.21037/atm-24-40/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://atm.amegroups.com/article/view/10.21037/atm-24-40/coif). Z.J.W. reports that stipend funding was provided by NIH/NCATS Colorado CTSA T32TR004366. L.M.P. reports that this study was funded by the Foundation for the Horse, Grayson Jockey Club Research Foundation, Colorado State University College Research Council Interdisciplinary Pilot Award, NIH/NCATS CTSA 5TL1TR002533-02 (L.M.P.), NIH 5T32OD010437-19 (L.M.P.), and Carolyn Quan and Porter Bennett. The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was approved by the Institutional Animal Care and Use Committee at Colorado State University (IACUC protocol #926) and conducted according to the national guidelines under which the institution operates and NIH Guidelines for the Care and Use of Laboratory Animals (8^th edition).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Cui A, Li H, Wang D, et al. Global, regional prevalence, incidence and risk factors of knee osteoarthritis in population-based studies. EClinicalMedicine 2020;29-30:100587. [Crossref] [PubMed]
2. Yelin E, Weinstein S, King T. The burden of musculoskeletal diseases in the United States. Semin Arthritis Rheum 2016;46:259-60. [Crossref] [PubMed]
3. Murphy L, Schwartz TA, Helmick CG, et al. Lifetime risk of symptomatic knee osteoarthritis. Arthritis Rheum 2008;59:1207-13. [Crossref] [PubMed]
4. Ireland JL, Clegg PD, McGowan CM, et al. Disease prevalence in geriatric horses in the United Kingdom: veterinary clinical assessment of 200 cases. Equine Vet J 2012;44:101-6. [Crossref] [PubMed]
5. Neundorf RH, Lowerison MB, Cruz AM, et al. Determination of the prevalence and severity of metacarpophalangeal joint osteoarthritis in Thoroughbred racehorses via quantitative macroscopic evaluation. Am J Vet Res 2010;71:1284-93. [Crossref] [PubMed]
6. Frisbie DD, Cross MW, McIlwraith CW. A comparative study of articular cartilage thickness in the stifle of animal species used in human pre-clinical studies compared to articular cartilage thickness in the human knee. Vet Comp Orthop Traumatol 2006;19:142-6. [Crossref] [PubMed]
7. McIlwraith CW, Fortier LA, Frisbie DD, et al. Equine Models of Articular Cartilage Repair. Cartilage 2011;2:317-26. [Crossref] [PubMed]
8. Chu CR, Szczodry M, Bruno S. Animal models for cartilage regeneration and repair. Tissue Eng Part B Rev 2010;16:105-15. [Crossref] [PubMed]
9. Reesink HL, Watts AE, Mohammed HO, et al. Lubricin/proteoglycan 4 increases in both experimental and naturally occurring equine osteoarthritis. Osteoarthritis Cartilage 2017;25:128-37. [Crossref] [PubMed]
10. Kramer WC, Hendricks KJ, Wang J. Pathogenetic mechanisms of posttraumatic osteoarthritis: opportunities for early intervention. Int J Clin Exp Med 2011;4:285-98. [PubMed]
11. Manferdini C, Paolella F, Gabusi E, et al. From osteoarthritic synovium to synovial-derived cells characterization: synovial macrophages are key effector cells. Arthritis Res Ther 2016;18:83. [Crossref] [PubMed]
12. Fahy N, de Vries-van Melle ML, Lehmann J, et al. Human osteoarthritic synovium impacts chondrogenic differentiation of mesenchymal stem cells via macrophage polarisation state. Osteoarthritis Cartilage 2014;22:1167-75. [Crossref] [PubMed]
13. Fichadiya A, Bertram KL, Ren G, et al. Characterizing heterogeneity in the response of synovial mesenchymal progenitor cells to synovial macrophages in normal individuals and patients with osteoarthritis. J Inflamm (Lond) 2016;13:12. [Crossref] [PubMed]
14. Bondeson J, Wainwright SD, Lauder S, et al. The role of synovial macrophages and macrophage-produced cytokines in driving aggrecanases, matrix metalloproteinases, and other destructive and inflammatory responses in osteoarthritis. Arthritis Res Ther 2006;8:R187. [Crossref] [PubMed]
15. Van Lent PL, Van den Hoek AE, Van den Bersselaar LA, et al. In vivo role of phagocytic synovial lining cells in onset of experimental arthritis. Am J Pathol 1993;143:1226-37. [PubMed]
16. Li YS, Luo W, Zhu SA, et al. T Cells in Osteoarthritis: Alterations and Beyond. Front Immunol 2017;8:356. [Crossref] [PubMed]
17. Keller LE, Tait Wojno ED, Begum L, et al. T Helper 17-Like Regulatory T Cells in Equine Synovial Fluid Are Associated With Disease Severity of Naturally Occurring Posttraumatic Osteoarthritis. Am J Sports Med 2023;51:1047-58. [Crossref] [PubMed]
18. Menarim BC, Gillis KH, Oliver A, et al. Autologous bone marrow mononuclear cells modulate joint homeostasis in an equine in vivo model of synovitis. FASEB J 2019;33:14337-53. [Crossref] [PubMed]
19. Menarim BC, Gillis KH, Oliver A, et al. Inflamed synovial fluid induces a homeostatic response in bone marrow mononuclear cells in vitro: Implications for joint therapy. FASEB J 2020;34:4430-44. [Crossref] [PubMed]
20. Gu Y, Hu Y, Zhang H, et al. Single-cell RNA sequencing in osteoarthritis. Cell Prolif 2023;56:e13517. [Crossref] [PubMed]
21. Rai MF, Collins KH, Lang A, et al. Three decades of advancements in osteoarthritis research: insights from transcriptomic, proteomic, and metabolomic studies. Osteoarthritis Cartilage 2024;32:385-97. [Crossref] [PubMed]
22. Ji Q, Zheng Y, Zhang G, et al. Single-cell RNA-seq analysis reveals the progression of human osteoarthritis. Ann Rheum Dis 2019;78:100-10. [Crossref] [PubMed]
23. Tan Z, Chen R, Lin H, et al. The Identification of Immune-Related Biomarkers for Osteoarthritis Immunotherapy Based on Single-Cell RNA Sequencing Analysis. Genet Res (Camb) 2023;2023:5574636. [Crossref] [PubMed]
24. Sebastian A, Hum NR, McCool JL, et al. Single-cell RNA-Seq reveals changes in immune landscape in post-traumatic osteoarthritis. Front Immunol 2022;13:938075. [Crossref] [PubMed]
25. Liu W, Chen Y, Zeng G, et al. Single-Cell Profiles of Age-Related Osteoarthritis Uncover Underlying Heterogeneity Associated With Disease Progression. Front Mol Biosci 2021;8:748360. [Crossref] [PubMed]
26. Mizoguchi F, Slowikowski K, Wei K, et al. Functionally distinct disease-associated fibroblast subsets in rheumatoid arthritis. Nat Commun 2018;9:789. [Crossref] [PubMed]
27. Zhang F, Wei K, Slowikowski K, et al. Defining inflammatory cell states in rheumatoid arthritis joint synovial tissues by integrating single-cell transcriptomics and mass cytometry. Nat Immunol 2019;20:928-42. [Crossref] [PubMed]
28. Lewis MJ, Barnes MR, Blighe K, et al. Molecular Portraits of Early Rheumatoid Arthritis Identify Clinical and Treatment Response Phenotypes. Cell Rep 2019;28:2455-2470.e5. [Crossref] [PubMed]
29. McIlwraith CW, Frisbie DD, Kawcak CE, et al. The OARSI histopathology initiative - recommendations for histological assessments of osteoarthritis in the horse. Osteoarthritis Cartilage 2010;18:S93-105. [Crossref] [PubMed]
30. Baxter GM, Stashak TS, Keegan KG. Examination for Lameness: History, Visual Exam, and Conformation. In: Baxter GM. editor. Adams and Stashak’s Lameness in Horses. John Wiley & Sons, Inc.; 2020:67-188.
31. Gillespie CC, Adams SB, Moore GE. Methods and Variables Associated with the Risk of Septic Arthritis Following Intra-Articular Injections in Horses: A Survey of Veterinarians. Vet Surg 2016;45:1071-6. [Crossref] [PubMed]
32. Hao Y, Hao S, Andersen-Nissen E, et al. Integrated analysis of multimodal single-cell data. Cell 2021;184:3573-3587.e29. [Crossref] [PubMed]
33. McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: Doublet Detection in Single-Cell RNA Sequencing Data Using Artificial Nearest Neighbors. Cell Syst 2019;8:329-337.e4. [Crossref] [PubMed]
34. Hafemeister C, Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol 2019;20:296. [Crossref] [PubMed]
35. Zappia L, Oshlack A. Clustering trees: a visualization for evaluating clusterings at multiple resolutions. Gigascience 2018;7:giy083. [Crossref] [PubMed]
36. Aran D, Looney AP, Liu L, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol 2019;20:163-72. [Crossref] [PubMed]
37. Patel RS, Tomlinson JE, Divers TJ, et al. Single-cell resolution landscape of equine peripheral blood mononuclear cells reveals diverse cell types including T-bet(+) B cells. BMC Biol 2021;19:13. [Crossref] [PubMed]
38. Miller SA, Policastro RA, Sriramkumar S, et al. LSD1 and Aberrant DNA Methylation Mediate Persistence of Enteroendocrine Progenitors That Support BRAF-Mutant Colorectal Cancer. Cancer Res 2021;81:3791-805. [Crossref] [PubMed]
39. Zheng L, Qin S, Si W, et al. Pan-cancer single-cell landscape of tumor-infiltrating T cells. Science 2021;374:abe6474. [Crossref] [PubMed]
40. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014;15:550. [Crossref] [PubMed]
41. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2021;2:100141. [Crossref] [PubMed]
42. Liberzon A, Subramanian A, Pinchback R, et al. Molecular signatures database (MSigDB) 3.0. Bioinformatics 2011;27:1739-40. [Crossref] [PubMed]
43. Jin S, Guerrero-Juarez CF, Zhang L, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun 2021;12:1088. [Crossref] [PubMed]
44. Jablonski KA, Amici SA, Webb LM, et al. Novel Markers to Delineate Murine M1 and M2 Macrophages. PLoS One 2015;10:e0145342. [Crossref] [PubMed]
45. Zhu L, Zhou X, Gu M, et al. Dapl1 controls NFATc2 activation to regulate CD8(+) T cell exhaustion and responses in chronic infection and cancer. Nat Cell Biol 2022;24:1165-76. [Crossref] [PubMed]
46. Riihimäki M, Fegraeus K, Nordlund J, et al. Single-cell transcriptomics delineates the immune cell landscape in equine lower airways and reveals upregulation of FKBP5 in horses with asthma. Sci Rep 2023;13:16261. [Crossref] [PubMed]
47. Zambrano-Zaragoza JF, Romo-Martínez EJ, Durán-Avelar Mde J, et al. Th17 cells in autoimmune and infectious diseases. Int J Inflam 2014;2014:651503. [Crossref] [PubMed]
48. Nishikoba N, Kumagai K, Kanmura S, et al. HGF-MET Signaling Shifts M1 Macrophages Toward an M2-Like Phenotype Through PI3K-Mediated Induction of Arginase-1 Expression. Front Immunol 2020;11:2135. [Crossref] [PubMed]
49. Du L, Lin L, Li Q, et al. IGF-2 Preprograms Maturing Macrophages to Acquire Oxidative Phosphorylation-Dependent Anti-inflammatory Properties. Cell Metab 2019;29:1363-1375.e8. [Crossref] [PubMed]
50. Waring P, Müllbacher A. Cell death induced by the Fas/Fas ligand pathway and its role in pathology. Immunol Cell Biol 1999;77:312-7. [Crossref] [PubMed]
51. Curat CA, Wegner V, Sengenès C, et al. Macrophages in human visceral adipose tissue: increased accumulation in obesity and a source of resistin and visfatin. Diabetologia 2006;49:744-7. [Crossref] [PubMed]
52. Yun MR, Seo JM, Park HY. Visfatin contributes to the differentiation of monocytes into macrophages through the differential regulation of inflammatory cytokines in THP-1 cells. Cell Signal 2014;26:705-15. [Crossref] [PubMed]
53. Ni R, Guo XE, Yan C, et al. Hemodynamic stress shapes subchondral bone in osteoarthritis: An emerging hypothesis. J Orthop Translat 2022;32:85-90. [Crossref] [PubMed]
54. Findlay DM. Vascular pathology and osteoarthritis. Rheumatology (Oxford) 2007;46:1763-8. [Crossref] [PubMed]
55. McIlwraith CW, Frisbie DD, Kawcak CE. The horse as a model of naturally occurring osteoarthritis. Bone Joint Res 2012;1:297-309. [Crossref] [PubMed]
56. Sage SE, Leeb T, Jagannathan V, et al. Single-cell profiling of bronchoalveolar cells reveals a Th17 signature in neutrophilic severe equine asthma. Immunology 2024;171:549-65. [Crossref] [PubMed]
57. Gordon S, Plüddemann A, Martinez Estrada F. Macrophage heterogeneity in tissues: phenotypic diversity and functions. Immunol Rev 2014;262:36-55. [Crossref] [PubMed]
58. Croft AP, Campos J, Jansen K, et al. Distinct fibroblast subsets drive inflammation and damage in arthritis. Nature 2019;570:246-51. [Crossref] [PubMed]
59. Liu B, Zhang M, Zhao J, et al. Imbalance of M1/M2 macrophages is linked to severity level of knee osteoarthritis. Exp Ther Med 2018;16:5009-14. [Crossref] [PubMed]
60. Collin M, Bigley V. Human dendritic cell subsets: an update. Immunology 2018;154:3-20. [Crossref] [PubMed]
61. Xiao J, Zhang P, Cai FL, et al. IL-17 in osteoarthritis: A narrative review. Open Life Sci 2023;18:20220747. [Crossref] [PubMed]
62. Liang Y, Lin F, Huang Y. Identification of Biomarkers Associated with Diagnosis of Osteoarthritis Patients Based on Bioinformatics and Machine Learning. J Immunol Res 2022;2022:5600190. [Crossref] [PubMed]
63. Faust HJ, Zhang H, Han J, et al. IL-17 and immunologically induced senescence regulate response to injury in osteoarthritis. J Clin Invest 2020;130:5493-507. [Crossref] [PubMed]
64. Liang D, Zuo A, Shao H, et al. IL-23 receptor expression on γδ T cells correlates with their enhancing or suppressive effects on autoreactive T cells in experimental autoimmune uveitis. J Immunol 2013;191:1118-25. [Crossref] [PubMed]
65. Lubberts E, Koenders MI, Oppers-Walgreen B, et al. Treatment with a neutralizing anti-murine interleukin-17 antibody after the onset of collagen-induced arthritis reduces joint inflammation, cartilage destruction, and bone erosion. Arthritis Rheum 2004;50:650-9. [Crossref] [PubMed]
66. Kemble S, Croft AP. Critical Role of Synovial Tissue-Resident Macrophage and Fibroblast Subsets in the Persistence of Joint Inflammation. Front Immunol 2021;12:715894. [Crossref] [PubMed]
67. Culemann S, Grüneboom A, Nicolás-Ávila JÁ, et al. Locally renewing resident synovial macrophages provide a protective barrier for the joint. Nature 2019;572:670-5. [Crossref] [PubMed]
68. Colbath AC, Dow SW, Hopkins LS, et al. Induction of Synovitis Using Interleukin-1 Beta: Are There Differences in the Response of Middle Carpal Joint Compared to the Tibiotarsal Joint? Front Vet Sci 2018;5:208. [Crossref] [PubMed]