
==== Front
iScience
iScience
iScience
2589-0042
Elsevier

S2589-0042(24)02052-2
10.1016/j.isci.2024.110827
110827
Article
Single-cell transcriptome and crosstalk analysis reveals immune alterations and key pathways in the bone marrow of knee OA patients
Chatterjee Paramita 12
Stevens Hazel Y. 12
Kippner Linda E. 12
Bowles-Welch Annie C. 12
Drissi Hicham 3
Mautner Kenneth 3
Yeago Carolyn 2
Gibson Greg greg.gibson@biology.gatech.edu
4∗
Roy Krishnendu krish.roy@vanderbilt.edu
15678∗∗
1 Marcus Center for Therapeutic Cell Characterization and Manufacturing, The Parker H. Petit Institute for Bioengineering and Biosciences, Georgia Institute of Technology, Atlanta, GA, USA
2 The Parker H. Petit Institute for Bioengineering and Biosciences Georgia Institute of Technology, Atlanta, GA, USA
3 Department of Orthopaedics, Emory University School of Medicine, Atlanta, GA 30322, USA
4 School of Biological Sciences, Georgia Institute of Technology, Atlanta, GA 30332, USA
5 Department of Biomedical Engineering, School of Engineering, Vanderbilt University, Nashville, TN, USA
6 Department of Pathology, Microbiology and Immunology, School of Medicine, Vanderbilt University, Nashville, TN, USA
7 Department of Chemical and Biomolecular Engineering, School of Engineering, Vanderbilt University, Nashville, TN, USA
∗ Corresponding author greg.gibson@biology.gatech.edu
∗∗ Corresponding author krish.roy@vanderbilt.edu
8 Lead contact

26 8 2024
20 9 2024
26 8 2024
27 9 11082714 11 2023
10 6 2024
22 8 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Summary

Knee osteoarthritis (OA) is a significant medical and economic burden. To understand systemic immune effects, we performed deep exploration of bone marrow aspirate concentrates (BMACs) from knee-OA patients via single-cell RNA sequencing and proteomic analyses from a randomized clinical trial (MILES: NCT03818737). We found significant cellular and immune alterations in the bone marrow, specifically in MSCs, T cells and NK cells, along with changes in intra-tissue cellular crosstalk during OA progression. Unlike previous studies focusing on injury sites or peripheral blood, our probe into the bone marrow—an inflammation and immune regulation hub—highlights remote organ impact of OA, identifying cell types and pathways for potential therapeutic targeting. Our findings highlight increased cellular senescence and inflammatory pathways, revealing key upstream genes, transcription factors, and ligands. Additionally, we identified significant enrichment in key biological pathways like PI3-AKT-mTOR signaling and IFN responses, showing their potentially crucial role in OA onset and progression.

Graphical abstract

Highlights

• OA affects immune homeostasis and has systemic effects on bone marrow

• Key interactions between immune and stromal cells in OA bone marrow is revealed

• MSCs in OA bone marrow show most gene changes and strong immune cell communication

• Pathway analysis points to altered autophagy and inflammation in OA

Biopsy sample; Computational bioinformatics; Transcriptomics

Subject areas

Biopsy sample
Computational bioinformatics
Transcriptomics
Published: August 26, 2024
==== Body
pmcIntroduction

Osteoarthritis (OA) is the most prevalent type of arthritis primarily affecting the knees, hips, hands, and spine, and is marked by joint degeneration.1,2,3,4,5 As the greatest global cause of pain and disability, it impacts over 32.5 million American adults alone, leading to an estimated yearly healthcare cost of $27 billion. Primary risk factors for OA are aging and obesity,6,7,8 which imply that excessive “wear and tear” and amplified biomechanical stress leads to joint degeneration. However, the relationship between OA and other health factors is complex and cannot be solely explained by me5chanical pressure, e.g., the higher occurrence of OA in the non-weight-bearing joints of obese individuals. To date, there are no licensed disease-modifying OA drugs, but non-surgical alternatives, such as orthobiologics, are attractive treatment candidates currently being evaluated in clinical trials.9,10,11,12

Emerging research points to a persistent mild inflammatory condition in OA, resulting in a systemic release of pro-inflammatory substances like adipokines from fat tissue, which may play a larger role in obesity-related OA.11 Changes in cellular structure in distant or seemingly unrelated bodily tissues caused by OA are being recognized as potential OA correlates. Twin studies further suggest that genetic factors account for more than half of the risk variance, and possibly up to 70% for spinal OA, with some of these factors likely having a systemic effect.12 Bone marrow aspirate concentrate (BMAC) is an autologous cell therapy typically harvested from multiple sites in the iliac crest. BMAC contains mesenchymal stromal cells (MSCs) that have been widely investigated for immunomodulatory and reparative effects.13 Typically, bone marrow aspirates contain very few MSCs (0.001–0.01%14) compared to the other cell populations; thus, procurement techniques to maximize the number of MSCs such as concentrating to produce BMAC are investigated as approaches to improve clinical outcomes.15 BMAC is also an enriched source of lymphocytes (13%), eosinophils (2.2%), monocytes (1.3%), basophils (0.1%),16 platelets (8.7-fold increase), and MSCs (0.03%) compared to bone marrow aspirate.17 Several cell-types whose abundance in circulation is correlated with OA pathology, for example, plasma cells, which originated from the bone marrow. Thus, a better understanding of the pathogenesis and disease severity of OA may be informed by deeper investigation into correlations between alterations identified in circulation and distal tissues, including the bone marrow.

The innate and adaptive immune systems play a key role in the inflammatory pathogenesis in OA. Evidence from the peripheral blood from individuals with OA has shown altered CD4+/CD8+ T cell ratios18,19 and pro-inflammatory T cell polarization in early knee OA, promoting helper T (Th) cell subsets, specifically Th1, Th9, and Th17.20 Similar pro-inflammatory polarization has been seen in macrophage populations, and correspondingly, NF-κB inhibitors have been shown to be effective against the progression of macrophage-mediated OA.21 Intermediate and non-classical monocytes derived from OA patients had elevated CD16+ and activation markers HLA-DR and CCR2.22 In a retrospective study of knee OA, decreased neutrophil/monocyte ratio in the blood correlated with disease severity (KL grade), highlighting the active recruitment of monocytes and macrophages to the synovium.23 Indeed, a microarray study of peripheral blood monocytes identified 791 upregulated and 440 downregulated differentially expressed genes (DEGs) in OA.24 Given the various cell types that reside in bone marrow, some of which are concentrated in BMAC, we have examined features of BMAC that may correlate with disease severity of patients with knee OA.

Understanding the bone marrow’s role in OA can be significantly advanced by leveraging modern single-cell technologies, which allow for unprecedented insights into cellular behaviors and interactions within this complex tissue. Single cell RNAseq (scRNAseq) has become a valuable tool to investigate high dimensional gene expression analysis of cells. In this study, we examined BMAC from patients with OA in the MILES (NCT03818737) trial which is one of the largest clinical trials investigating cell therapies for knee OA. We used scRNAseq to define the cellular composition of BMAC and identify gene expression attributes that are likely to contribute to their mode of action in OA, particularly in easing joint pain, and compared them with publicly available NIH study of scRNAseq data from a cohort of non-OA donors.25 In this study, we compare BMAC obtained via different isolation protocols, specifically the EmCyte system and Ficoll density gradient centrifugation. While these methods have inherent differences that could influence the cellular profiles observed, we have been careful to match phenotypically similar cell types across OA and non-OA datasets to ensure that our findings focus on biological rather than methodological differences. Despite these precautions, the differences in isolation techniques may introduce subtle biases in cell subset recovery and viability, which could affect the interpretation of our results. To enhance the reliability of our results, we ensured that the non-OA and OA samples were collected using meticulous batch correction, and utilized a separate OA dataset of 56 patient samples from the same trial cohort for validation of OA relevant results, confirming that identified pathways and signatures are biological signatures and not technical anomalies. We employed robust analytical pipelines to evaluate changes in cell population abundances, identify differentially expressed genes, and cell-cell signaling mechanisms in OA versus non-OA BMAC. The validation of key findings in the replication dataset as well as steps taken to make like-for-like comparisons of gene expression in the same cell types, highlights gene expression differences in the bone marrow environment, offering unique insights into osteoarthritis pathology in patient bone marrow.

Results

Utilizing samples collected during a randomized, double-blinded multicenter clinical trial (MILES: NCT03818737), high dimensional transcriptomic analysis and cellular phenotyping of BMAC at the single cell level from OA patients (n = 77) was performed. For discovery, 21 of these OA samples were randomly chosen to be quantitatively compared to a similarly sized non-OA cohort (n = 20) extracted from a public database.25 We used the remaining 56 samples as a validation dataset for reproducibility of the OA relevant signatures found from the comparison of non-OA and OA cohorts. The patient and generated data information is listed in the Table S1. Following the workflow schematized in Figure 1, we compared cell-type proportions and differential expression in bone marrow aspirate concentrate (BMAC) cells from non-OA and OA patients. Using scRNAseq data generated on the 10X Genomics Chromium platform followed by Illumina sequencing, we had 136,484 cells available for analysis after quality control. Henceforth, we refer to the healthy as non-OA and diseased as OA throughout the article for the convenience of the reader. After QC filtering to remove low quality cells and mitochondrial genes, and log normalization, the two datasets were integrated using CCA (canonical correlation analysis) based on disease status and then cluster generation based on the transcriptomic profile of the cells was performed to classify BMAC cell types (see STAR Methods for details).Figure 1 Graphical overview of single-cell omics analysis of patient OA and non-OA bone marrow

1–9 steps illustrating the experimental and analytical schematic process. The cryopreserved cells were processed to generate single cell RNA (scRNAseq) library generation of OA patient bone marrow samples. These samples were parallelly processed for mass cytometry and flow cytometry assessments. After sequencing the scRNAseq libraries, the cell types were identified from cluster generation. DEG, cellular crosstalk, and pathway analyses were performed to correlate and interpret the relevant biological findings.

Cell proportion and abundance of phenotypically similar cell types reveal increased abundance of NK cells, T cells, and megakaryocytes in OA

We next identified phenotypically similar cell types from both OA and non-OA datasets, facilitating accurate comparison of their abundances and proportions. After clustering, we identified 16 major cell types that were then used to evaluate the heterogeneity of the BMAC as a whole and, to determine which cell types may differ between the OA and non-OA groups. We further investigated the robustness of our comparative analysis by including an assessment of excluded erythrocytes in BMAC samples, which helped address potential discrepancies from differential cell exclusion during processing. Additionally, we analyzed both fresh and cryopreserved BMAC samples from the same donors, finding minimal variation in gene expression, which supports the robustness of our findings (Figure S11). DEG analysis was performed using the Wilcoxon rank-sum test with MAST to control for batch effects. The compiled analysis of all donor cells was annotated and then divided into two to show the OA and non-OA cell clusters side by side (Figure 2A).Figure 2 Single-cell transcriptomics profiling of BMAC cells

(A) Cell type identification of OA and non-OA BMAC cells with SEURAT version 4 analytical pipeline. The UMAP shows the cell populations for both groups split showing the cell type clusters. Both OA and nonOA groups have all the major cell types present.

(B) A stacked barplot shows the cell type cluster proportions comparing the diseased group (OA) and the non-diseased group (non-OA). The significant changes in the cell proportions are shown with the p value significance. The significance of p value is listed as '∗∗' for less than or equal to 0.01 and '∗∗∗' for less than or equal to 0.001. Overall, we observe an increase of NK cells, megakaryocytes, double-negative T cells (DNT), CD8 T cells, MSCs, and progenitor B cells in the OA cohort.

(C) Donor-wise distribution on cell cluster proportion to highlight the cell type variabilities in BMAC donors.

(D) Dotplot shows the expression of the key markers to identify cell type for each cell type. Red is nonOA and Blue is OA. The size of the dots represents the expression percentages.

Analysis of cell proportions revealed that NK cells and megakaryocytes (MEGA) were significantly increased in the OA group (both p < 0.001) relative to the non-OA group, while double-negative T cells (DNT, p < 0.01), as well as CD8 T cells, MSCs, and B cells also showed this tendency. Conversely, there was depletion of CD14 and CD16 monocytes, progenitor dendritic cells, and plasma B cells (all p ≤ 0.001). These proportions are shown across all samples in Figure 2B, and by individual in Figure 2C, illustrating the heterogeneity within and between OA and non-OA. Details of the statistical analysis are included in supplementary document Table S2.

Mass cytometry analysis, summarized by representative viSNE projections in Figure 3A and Table S4, also revealed a significantly reduced total T cell population (p < 0.001) in OA BMAC compared to non-OA BMAC. Transcriptomes of OA BMAC show a significant decrease in CD4 T cells (p < 0.05) compared to non-OA BMAC, which is consistent with the marked reduction in CD4 T cells measured by mass (p < 001) and flow (p < 0.05) cytometry compared to non-OA BMAC. A complementary reduction in CD8 T cell frequency in OA BMAC was measured by both mass (p < 0.01) and flow (p < 0.001) cytometry analyses, however, this was not apparent in the scRNAseq analysis. Since CyToF and scRNAseq identify different cell types, pairwise comparison is not supported, yet similar overall trends are documented in Figure 3B. These include directional shifts in non-classical (CD16+) monocytes, DNT, and HSPC.Figure 3 Mass cytometry cell type identifications and abundance analysis of the OA BMAC samples

(A) Representative mass cytometry profile of the OA donors with the cell type annotated.

(B) The bar plots show cell-type abundances correlation in scRNAseq and mass cytometry side by side for both OA and nonOA cohorts. The error bars on each data point represent the standard error of the mean (SEM) for the cell type frequencies across biological replicates. The SEM provides an estimate of the variability in cell type abundances within each sample group.

Variance analysis reveals greater impact of disease status than BMI on gene expression across cell types in osteoarthritis

Since inflammation or other physiological responses attributable to high body weight might contribute to cellular function, we analyzed the variance components of gene expression attributable to BMI and disease status across different cell types (Figure S10). Since BMI measurement was not available in the non-OA dataset, this comparison was restricted to the OA dataset, but nevertheless allows benchmarking against the disease contrast. The variance partitioning results show that while BMI does contribute to variance in certain cell types such as MSCs and mono CD14 cells, the magnitude of the effect is much smaller than that due to disease status across all cell types (Figure S10A) and in most individual cell types, including CD8 T cells (Figure S10B), NK cells, and plasma B cells. The only exception, where BMI has a more noticeable impact, was the MSCs, consistent with body weight indirectly modulating the transcriptome and likely function of this potential mediator of therapeutic response (Figure S10C).

Cell communication analysis and crosstalk activity reveal enriched ligand-receptor pairs related to bone modulation and degeneration pathways in OA

We next used cell crosstalk analysis to determine which cell types are likely to drive differences between OA and non-OA datasets. Ribbons linking cell types in the plots in Figure 4A indicate inferred inter-cellular communication on the basis of the expression of receptors and ligands in the respective cell types. Comparison of OA and non-OA revealed several cell types with increased probable crosstalk, notably involving MSC, which may be signaling to 6 other cell types in non-OA but nine other cell types in OA. This result was replicated in the validation dataset. On the other hand, progenitor B cell interactions were not observed in the OA datasets. The top enriched ligand receptors for both OA and non-OA cohorts are highlighted with wordcloud visualization in Figure 4B. We also observed 18 ligand-receptor pathway enrichments that were exclusive to the OA cohort (blue text on the Y axis) and just 3 exclusive pathways significantly enriched (red text on the Y axis) in non-OA (Figure 4C). The list of enriched ligand-receptors is provided in the supplementary (Table S3).Figure 4 Detailed cellular crosstalk analysis in BMAC samples

(A) Cellular interaction profiling using CellChat.

This panel visualizes the top 20 cellular interactions identified within the OA cohort, across both training and validation datasets. The analysis highlights dynamic communication networks, predominantly featuring interactions between mesenchymal stem cells (MSCs) and other critical immune cells such as NK cells, T cells, and hematopoietic stem and progenitor cells (HSPCs). The graphical representation delineates the complexity of intercellular communications, underlining the enriched signaling pathways that potentially influence therapeutic outcomes.

(B) The word cloud illustrates the relative prominence of ligand-receptor pairs in the OA and non-OA groups, emphasizing the differential expression of key molecules. Notable ligands such as LAMA4, LAMB1, BMP5, LAMC1, PTPRM, and SEMA4A are exclusively enriched in the OA cohort, suggesting a unique molecular signature that may be pivotal in OA pathogenesis and progression.

(C) Comparative analysis of enriched ligand pathways: presented as a bar plot, this panel quantifies and compares the enriched ligand-receptor pathways between OA and non-OA groups. Color-coded for intuitive interpretation (OA in Blue, non-OA in Red), the plot provides a visual summary of the pathway distribution, highlighting the presence of multiple ligand-receptor interactions unique to the OA group. This differential pathway activity could inform targeted therapeutic strategies.

From the differential crosstalk analysis,26 we further identified the top ligand-receptor pairs which may be engaged in OA-BMAC intercellular signaling and warrant further investigation. Notably, the laminins, LAMC1, LAMA4, LAMB1, as well as growth factor receptors FN1, BMP5, and SEMA4D are among the top enriched ligand receptors in OA (Figures 4B and 4C). These three laminins have been shown to influence hypertrophic chondrocyte clustering in OA development.27 Also, the fibronectin 1 protein encoded by FN1 has been found to stimulate the production of matrix metalloproteinases (MMPs) in cartilage which can contribute to the degradation of cartilage in OA.28 Similar ligand receptor pairs were enriched in the validation OA dataset, strengthening our analytical findings related to observed correlates of OA pathogenesis.

Interestingly, SEMA4D emerged as one of the ligands showing high activity in MSCs in OA (Figure 5A). Two previous studies have reported the expression of SEMA4D in RANKL-activated osteoclasts but not in osteoblasts.29,30 Along with SEMA4D we also found that the cytokine IL-16 is highly enriched in OA, in both the training and validation datasets (Figure 5A). It plays a role in the inflammatory response by recruiting and activating immune cells through activation of MMP genes.31,32Figure 5 Analysis of key pathways and gene expressions enriched in OA Using BMAC from OA patients

(A) This panel illustrates the ligand-receptor activity specifically associated with mesenchymal stem cells (MSCs) under the OA condition, utilizing Nichenet analysis to pinpoint critical interactions. It shows a detailed comparison across non-OA, OA training, and OA validation datasets, highlighting several key ligand-receptor pairs. Notable among these are antigen-presenting receptors and SEMA4D, which are consistently enriched in both OA datasets, indicating their potential role in the molecular mechanisms underpinning OA.

(B) The bar plots in this panel depict the intensity and distribution of IL16 and SEMA4 ligand pathways, revealing significant communication patterns within the MSCs across different study groups. This comparative analysis underscores the heightened activity in both the OA training and validation datasets relative to non-OA controls. The graphical representation provides insights into the differential expression and potential functional implications of these pathways in OA pathophysiology.

Given the enriched ligand-receptor information, we then used Nichenet to further investigate the possible target genes for these selected ligand-receptor pathways33 and subsequently DEGs for the downstream pathway analysis (Figure S1). Differential cell abundances were also confirmed using the Nichenet tool (Figure S2). From Overall communication analysis, we found that ligands and their target genes participating in the RANK/RANKL/OPG,34 cGAS-STING, and cellular senescence pathways are highly expressed in the OA dataset. It is possible that the catabolism observed in human OA cartilage, which is mediated by MSCs obtained from patients with OA, could be attributed to TGF-β lateral signaling activating BMP.35,36 MSCs appear to communicate most extensively with NK cells, naive T cells, and HSPCs. We also found that there are several ligand-receptor enrichments present only in OA along with the ligands significantly more highly enriched in OA, as summarized for IL16 and SEMA4 and MHC (HLA) in Figures 5A and 5B. A large portion of these ligands including the MHC-I also showed interactions with the MSC populations (Figure 5A).

Differential gene expression in BMAC cell types indicate highest number of differentially expressed genes in the OA MSC population

We next contrasted the gene expression in BMAC cell types from both groups, delving deeper into the cell types showing increased communication within the osteoarthritis (OA) group. Using further subset analysis, we intended to identify and understand the biologically pertinent characteristics associated with OA. We identify the population as MSC which has CD90 (THY1), CD73 (NT5E), and CD105 (ENG) positive expression while they lack CD45 (PTPRC), and CD34 expression (Figure S12). These are the MSC markers as per ISCT guidelines (Dominici, 2006 #198).

Interestingly, the greatest number of differentially expressed genes was observed in the MSC population (2082 genes) (Figure 6A). Expression of key signaling molecules in MSC is summarized in Figure 6B, including four genes, CDKN1A, CYR61, IL10RB, and IL7, whose expression is highly elevated in the non-OA group and minimally detected in the OA group (Figure 6B). CDKN1A deficiency is related to susceptibility to enhanced inflammation in OA,37 which pushes the immune system toward the induced NF-κB pathway. CYR61 regulates cell adhesion and migration and is required to maintain the properties of bone marrow MSCs.38 The decrease of CYR61 in OA suggests that retention of proliferation and growth factor responsiveness is altered in this group. One study found RA pathogenesis promoted via the IL17 pathway39 but this remains to be investigated with respect to OA pathogenesis. Our data suggest that CYR61 downregulation poses a vital role in chronic inflammation in OA pathogenesis and can be a potential therapeutic target.Figure 6 Comprehensive analysis of differential gene expression in OA

(A) This panel highlights the results of differential expression analysis between OA and non-OA groups, specifically within mesenchymal stem cells (MSCs). It showcases a bar graph depicting the relative number of differentially expressed genes (DEGs), where MSCs exhibit the highest disparity in gene expression, underscoring their pivotal role in OA pathophysiology.

(B) This section delves into the specific genes within MSCs that are altered in OA compared to non-OA conditions. The analysis identifies key genes implicated in the pathogenesis of OA, offering insights into the molecular alterations that may drive disease progression. This visualization aids in understanding the gene-level changes and their potential impact on therapeutic targeting.

(C) Focusing on toll-like receptors, this panel presents a clear visualization of their differential expression levels in the OA dataset relative to non-OA. The significant upregulation of these receptors in OA suggests their crucial involvement in inflammatory responses and innate immunity within the OA context, potentially contributing to disease mechanisms and symptoms.

Analysis of the T cell subsets revealed that CD8 T cells and naive T cells were highly activated in the OA compared to the non-OA group. Pseudotime analysis suggests that T cells differentiate at an earlier time point than B cells in OA, whereas T cells were the final trajectory point in non-OA samples (Figures S3A and S3B). This analysis suggests that there may be early recruitment of T cells in OA. Additionally, IL17RA expression was particularly high in the OA compared to non-OA in various T cell subsets (Figure S4), which is notable since previous studies have shown that IL17RA contributes to OA pathophysiology.40

Multiple other cell types display interesting differential gene expression. NK cells, which had increased abundance in the OA group and communicated directly with MSCs, were subdivided into 16 sub-clusters implicating different states of NK activation between OA and non-OA (Figure S5). Classes of differentially expressed genes that are likely engaged by MSC polarization and licensing include a Toll-like receptor cascade, cellular senescence, and osteoclast differentiation pathways, also implicating TLR4 in OA cartilage degradation. Megakaryocytes were significantly increased in the OA group, with high expression of PF4 (encoding the cytokine platelet factor 4) and multiple genes associated with TLR activation and arthritis pathways (Figure S6). Mature B cells showed differentially expressed genes involved in oxidative stress, TLR regulation, and response to lipopolysaccharides (Figure S7). These pathways are associated with synovial inflammation, cartilage degradation, and matrix synthesis. DNT cells showed increased expression of LGALS3 and genes involved in cartilage degradation, chondrocyte homeostasis, and cellular senescence (Figure S8). Finally, monocytes exhibited gene expression changes indicative of cellular senescence. CD16+ monocytes, which were observed to be decreased in OA overall, had upregulated gene expression related to apoptosis, chemokines, focal adhesion, metabolic pathways, and osteoclast differentiation. The OA cohort showed elevated gene expression for transcription factors involved in maturation and activation.

OA-relevant signaling pathways identified through gene set enrichment and pathway analysis highlight changes in autophagy, inflammation, and senescence

In an attempt to elucidate the central pathways implicated in osteoarthritis pathogenesis and prognosis, we conducted a comprehensive gene ontology analysis employing the gene set enrichment analysis (GSEA) tool along with the Molecular Signatures Database (MSIGDB) hallmark gene sets. This analysis revealed several key pathways, including P53, PI3/AKT/MTOR signaling, oxidative phosphorylation, TGF-Beta signaling, IL2-STAT5 signaling, and interferon alpha and gamma responses, that showed substantial enrichment in OA (Figure 7A). Involvement of the unfolded protein response pathway is notable since endoplasmic reticulum stress may be a feature of osteoarthritic chondrocytes (Figure S9). Activation of this pathway implicates inflammation and apoptosis, potentially instigating cartilage degeneration.Figure 7 Pathway enrichment analysis

(A). GSEA pathways enriched in OA from the 20-sample cohort. Positive NES in red and negative NES in blue. The X axis shows the enriched pathways in OA dataset and Y axis shows the cell types. The heatmap color intensity signifies the NES scores for each pathway and celltypes in OA.

(B) GSEA pathways enriched in OA from the 56-sample validation cohort. Positive NES in red and negative NES in blue. The X axis shows the enriched pathways in OA dataset and Y axis shows the cell types. The heatmap color intensity signifies the NES scores for each pathway and cell types in OA.

Further investigations suggest that these identified pathways are crucial for the regulation of autophagy and cellular senescence, which are fundamental facets of the immune response associated with OA. Importantly, autophagy is the inherent process of eliminating senescent cells to maintaining tissue homeostasis. Disruption of autophagy was observed to instigate inflammation driven by senescent cells in OA.41 Moreover, we found that autophagy significantly contributes to the regenerative abilities of MSC, and therefore, its inhibition could negatively affect these cell properties. MSCs possibly playing a key role in stimulating further recruitment of immune cells, thereby exacerbating inflammatory pathways, cellular senescence, and inflammation. Contrastingly, such a characteristic pro-inflammatory signature of MSCs was not observed in the non-OA cohort.

Additionally, overview of the GSEA analysis showed substantial enrichment of the p53, PI3/AKT/mTOR, IL6/JAK/STAT, mTORC, and interferon signaling pathways in OA.39,42,43,44,45 Similarly, perturbation pathway scoring with the r package Progeny46 contrasting the OA and non-OA overall datasets identified estrogen, TGFB, hypoxia, and P53 pathways. We also found that the WNT pathway is very much downregulated in OA compared to non-OA datasets, and a similar signature was observed in the validation OA dataset as well (Figure 7B).

Discussion

This study presents a comprehensive view of the systemic knee-osteoarthritis (OA) pathology at the cellular level, utilizing single-cell transcriptomics to compare OA and healthy bone marrow. Previous studies related to omics landscape in OA has focused on OA relevant tissues like blood, plasma, urine, local joint cartilage, and synovium.47 Our methodological approach underlines the value of single-cell genomics in the field of rheumatology, providing a deeper understanding of the cellular diversity and functional states in OA.

In this study, we observed altered expression profiles and communication signatures in key cell types including NK cells, MSCs, monocytes, B cells, and T cells. Notably, MSCs exhibited a distinctly different signature in OA data. This difference suggests the presence of immunomodulatory changes within the patient’s bone marrow cells. It is plausible that these alterations could be the driving force behind the inflammatory signatures observed in knee-OA. One caveat to this study is that our control non-OA dataset is curated from a NIH public data source. To mitigate this limitation, we have taken four steps. First, we ensured that the BMAC isolation method was performed following a similar protocol with our OA datasets. Second, batch correction on the groups and conditions was performed, and preliminary data integration method comparison evaluated to find the most suitable and reproducible pipeline. This process revealed that even though our OA samples were collected from five different sites, the difference between the OA and non-OA conditions was much greater than among the sites. Third, we included a larger separate OA dataset of 56 patients as a validation dataset from our same OA cohort to replicate the key findings. Lastly, we observed that all of the identified pathways and signatures are biologically relevant and plausibly associated with knee-OA disease, an unlikely outcome for purely technical sources of error. These steps increase confidence that the study findings represent a true investigation of the perturbations that characterize bone marrow of osteoarthritis patients.

Previous studies have primarily focused on the injury site or peripheral blood mononuclear cells.48 Since our study uniquely investigates cells from the bone marrow, which is a systemic reservoir of inflammation and immune regulation,49 in a high-resolution and detailed manner, it highlights novel aspects of the systemic response to OA. Consistent with the OA literature, our data reveal increased cellular senescence and inflammatory pathways in OA, but we also identified dominant cell types and master pathways likely involved in the intercellular crosstalk in OA. These include key upstream genes that could serve as potential therapeutic targets for modulation or inhibition. With the growing interest in BMAC as a cell therapy for OA, our findings also suggest critical quality attributes (CQAs) for BMAC. The full therapeutic potential and implications of this approach warrant further exploration.

The observed T cell depletion in OA prompts several hypotheses, including senescence, exhaustion, or depletion. Future studies are encouraged to examine specific gene markers that define these states in our dataset, potentially leading to novel therapeutic targets. From previous studies50 senescent T cells communicate with non-lymphoid cells in pro-inflammatory conditions, and this can lead to tissue damage. Primary markers for senescence in CD8 T cells are CD57 (B3GAT1) and KLRG1 and loss of CD28 markers. Our analyses find high expression of the B3GAT1 and KLRG1 genes and loss of CD28 in OA cell groups, consistent with the inference of cellular senescence in the inflammatory conditions of the OA bone marrow. A combination of in vitro and in vivo experimentation has shown51 that Th17 cells induce senescence in fibroblasts and skew naive T cells toward Th17 or Th1 depending on the presence of TGFβ. By blocking IL17 activity, joint degradation was reduced and the marker CDKN1A was downregulated. In our study, we see a very high expression of CDKN1A in OA compared to the healthy non-OA group, consistently suggesting that IL17 and upstream mTORC1 pathway-associated markers can be used as signatures for OA.37 Previous study-based flow cytometry on blood from osteoarthritis patients has also found T cell and B cell subsets involved in the pathogenicity of OA.5

Our study utilized cell crosstalk analysis to pinpoint which cell types might drive differences between the OA and non-OA datasets. We observed notable increases in probable crosstalk involving mesenchymal stem cells (MSCs), which signal to a greater number of other cell types in the bone marrow in OA than in non-OA. This differential crosstalk highlights the potential role of MSCs in modulating the bone environment in OA, particularly through the activation of pathways related to bone modulation and degeneration. Although previously senescent MSCs has been reported in subchondral and chondrocyte studies,52 our study provides novel insight into the outgoing signaling ligand receptors which are involved in this crosstalk.

Key ligand-receptor pairs identified, such as laminins (LAMC1, LAMA4, LAMB1), fibronectin 1 (FN1), and SEMA4D, underscore the involvement of these molecules in OA pathophysiology. Laminins, for example, have been linked to hypertrophic chondrocyte clustering in OA development, suggesting a mechanism by which these extracellular matrix components could influence chondrocyte behavior and contribute to disease progression. FN1’s role in stimulating matrix metalloproteinases hints at its involvement in the degradation processes characteristic of OA. SEMA4D’s heightened activity in OA-specific MSCs points toward its role in immune modulation and bone remodeling, further supported by its expression in osteoclasts but not osteoblasts.

The differential expression and enriched activity of these ligand-receptor pairs, replicated in our validation dataset, lend strong support to our analytical findings and suggest that these molecules are not merely bystanders but active participants in OA pathogenesis. This is further supported by the analysis of target genes for these ligand-receptor pathways, revealing involvement in critical pathways such as RANK/RANKL/OPG, cGAS-STING, and cellular senescence, all of which are pivotal in OA’s inflammatory and degenerative processes.

Significant enrichment in several biological pathways in knee osteoarthritis was detected using gene set enrichment analysis (GSEA). Pathways such as PI3-AKT-mTOR signaling, IFN-gamma response, IFN-alpha response, unfolded protein response, TGFB signaling, IL2-STAT5 signaling, IL6-JAK-STAT3 signaling, and mTORC1 signaling, could potentially underpin key molecular mechanisms contributing to the onset and progression of osteoarthritis. The PI3-AKT-mTOR signaling pathway plays a crucial role in cell survival, proliferation, and differentiation. In the context of knee osteoarthritis, dysregulated PI3-AKT-mTOR signaling could influence chondrocyte survival and cartilage matrix synthesis, thereby contributing to disease pathogenesis.53 Interferon (IFN) responses, both IFN-gamma and IFN-alpha, are also enriched in our analysis. These pathways mediate immune responses and have been implicated in inflammatory diseases. As osteoarthritis has an inflammatory component, particularly in synovial inflammation, it is plausible that these pathways could play a role in the disease process. The unfolded protein response pathway is involved in endoplasmic reticulum stress response, a process shown to be increased in osteoarthritic chondrocytes. Its activation may lead to inflammation and apoptosis, driving cartilage degeneration. TGFB signaling plays a complex role in cartilage homeostasis and osteoarthritis. On one hand, TGFB is crucial for maintaining cartilage integrity; while on the other hand, dysregulated TGFB signaling can drive osteoarthritis progression by promoting abnormal chondrocyte differentiation and extracellular matrix degradation. The IL2-STAT5 signaling pathway has a critical role in T cell proliferation and differentiation,54 which are pivotal in immune response regulation. Although traditionally considered a non-inflammatory condition, recent studies indicate the involvement of immune pathways in osteoarthritis, suggesting a potential role for IL2-STAT555 signaling. The IL6-JAK-STAT3 signaling pathway mediates various cellular processes,42 including inflammation. IL6 is known to be upregulated in osteoarthritis, leading to increased inflammation and cartilage degradation. Therefore, this pathway likely contributes to the progression of knee osteoarthritis. Lastly, the mTORC1 signaling pathway, an integral part of the larger mTOR pathway, is vital in regulating cell growth and metabolism. Dysregulation of this pathway may affect chondrocyte homeostasis and thus, contribute to osteoarthritis pathogenesis.43

Another notable highlight of our bioinformatic analysis is the apparently exaggerated role of MSC in inflammation, which as discussed by Grandi et al.,56 has potential implications for MSC-based therapies in OA. From earlier studies, it is known that high levels of TGF-β upregulate the RUNX2 and SMAD pathways, ultimately leading to bone degeneration.57 Also shown in Figure 5B, the top differentially expressed genes in the OA MSCs are RUNX2, RACK1, TEAD1, CRH, BMPR2, HGF, MMP14, ROCK1, THBS1, FGFR1, TGFBR3, and MMP2. These genes are involved in the induction of osteoarthritis,58,59 stimulation of inflammation,60 cartilage homeostasis,61 cartilage loss, and bone resorption pathways.38,39,42,43,58,60,62,63,64,65,66,67,68 We also find IL7 is highly expressed in non-OA compared to OA overall and specifically in MSCs. A recent study showed that IL7 plays a vital role in CAR-T cell efficiency,10 and it is also important for diabetic wound healing.69 A possibility is that modified MSCs promote wound healing. With this functional in vitro evidence, IL7 can be a novel candidate to modulate the MSCs in OA toward more regenerative and repair pathways.69,70 We also see that MSCs in OA have very high expression of BMP pathways which are crucial for bone development. A mouse study71 showed that BMP5 silencing inhibits chondrocyte senescence and apoptosis as well as osteoarthritis progression. BMP5 regulates chondrocyte senescence via the p38/ERK signaling pathway and BMP silencing in chondrocytes reduces senescence-associated secretory phenotypes in the knee joint. BMPR1A expression also regulates the ability of MSCs to bring about adipogenesis.72 Therefore, BMP5 and associated genes in the signaling pathways may also be a potential target to explore.

Our findings, coupled with the enriched pathways observed, suggest a multifaceted mechanism of OA involving altered cell signaling and intercellular communication, primarily mediated by MSCs and their interaction with other immune cells like NK cells, HSPCs, and T cells. These interactions likely contribute to the inflammatory milieu and the subsequent tissue degradation seen in OA and may be modified by obesity given the effect of BMI status on the MSC transcriptome.

In conclusion, the specificity of these changes to OA, along with their validation across datasets and their biological plausibility, strengthens the hypothesis that targeting these pathways could offer new therapeutic avenues for OA treatment. Future studies should focus on delineating these interactions further and exploring their therapeutic potential in clinical settings.

Our research specifically addresses the global impact of osteoarthritis (OA) by highlighting cellular and molecular changes in bone marrow samples from OA patients. This highlights not only the well-understood local effects of OA on subchondral bone marrow, synovium, synovial fluid, and cartilage, but also underscores the unresolved issues in more distant tissues such as blood and iliac crest bone marrow. Thus, our findings offer novel and valuable insights into the broader systemic effects of OA.

Limitations of the study

In this study, we compared BMAC obtained through different isolation protocols, specifically the EmCyte system and Ficoll density gradient centrifugation. While these methods inherently differ, potentially influencing the observed cellular profiles, we carefully matched phenotypically similar cell types across osteoarthritis (OA) and non-OA datasets to ensure our findings emphasize biological rather than methodological differences. Despite these precautions, the isolation technique disparities might introduce subtle biases in cell subset recovery and viability, potentially affecting our cell type abundance result interpretations.

To address this limitation, the validation of key findings was performed in the validation and replication dataset, coupled with steps taken to ensure like-for-like comparisons of gene expression in the same cell types. This highlights gene expression differences within the bone marrow environment. These differences offer unique insights into osteoarthritis pathology in patient bone marrow.

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact Krish Roy (krish.roy@vanderbilt.edu).

Materials availability

This study did not generate new unique reagents.

Data and code availability

• Single cell RNA-sequencing data are deposited in GEO (GSE274018) and SRA (PRJNA1144164). These data are not publicly available at the time of publication due to a clinical trial data embargo until all results are reported. Until release, data are available by request to the lead contact and/or corresponding authors.

• Original code can be found at https://github.com/pchatterjee7/BMAC_OA_nonOA.

• All other items regarding this study can be made available by requesting from the lead contact.

Acknowledgments

We wish to acknowledge the Cellular Analysis and Cytometry Core and the Molecular Evolution Core facilities at the Parker H. 10.13039/100008267 Petit Institute for Bioengineering and Bioscience at the 10.13039/100006778 Georgia Institute of Technology for the use of their shared equipment, services, and expertise. Funding for this work was provided by The Billie and Bernie 10.13039/100019408 Marcus Foundation to 10.13039/100006939 Emory University and the 10.13039/100006778 Georgia Institute of Technology , through the Marcus Center for Therapeutic Cell Characterization and Manufacturing (MC3M). K.R. was supported by the Robert A. Milton Endowed Chair. MC3M is also supported by the 10.13039/100017015 Georgia Tech Foundation and the 10.13039/100008065 Georgia Research Alliance .

Author contributions

P.C. conceived of the study; P.C., L.K., H.Y.S., and A.B.W. performed the experiments; P.C., H.Y.S., A.B.W., L.K., and A.M. analyzed the data; P.C., H.Y.S., A.B.W., and L.K. wrote the manuscript; K.M. provided BMAC samples; P.C., H.Y.S., A.B.W., L.K., G.G., C.Y., K.R. guided, reviewed and edited manuscript drafts. All authors reviewed and approved the final manuscript.

Declaration of interests

The authors declare no competing interests.

STAR★Methods

Key resources table

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
Antibodies	
	
Anti-Human CD3 (UCHT1)-170Er	Fluidigm Corporation	Cat#3170001B	
Anti-Human CD45 (HI30)-89Y	Fluidigm Corporation	Cat#3089003B	
Anti-Human CD11c (Bu15)-147Sm	Fluidigm Corporation	Cat#3147008B	
Anti-Human CD34 (581)-166Er	Fluidigm Corporation	Cat#3166012B	
Anti-Human CD146 (P1H12)-155Gd	Fluidigm Corporation	Cat#3155006B	
Anti-Human CD56 (NCAM16.2)-176Yb	Fluidigm Corporation	Cat#3176008B	
Anti-Human CD31 (WM59)-144ND	Fluidigm Corporation	Cat#3144023B	
Anti-Human CD25 (2A3)-149Sm	Fluidigm Corporation	Cat#3149010B	
Anti-Human CD4 (RPA-T4)-145ND	Fluidigm Corporation	Cat#3145001B	
Anti-Human CD8 (RPA-T8)-146ND	Fluidigm Corporation	Cat#3146001B	
Anti-Human CD19 (HIB19)-142ND	Fluidigm Corporation	Cat#3142001B	
Anti-Human CD14 (M5E2)-175Lu	Fluidigm Corporation	Cat#3175015B	
Anti-Human CD16 (3G8)-148ND	Fluidigm Corporation	Cat#3148004B	
Anti-Human HLA- DR (L243)-173Yb	Fluidigm Corporation	Cat#3173005B	
Anti-Human CD11b/Mac-1 (ICRF44)-209Bi	Fluidigm Corporation	Cat#3209003B	
Anti-Human CD123/IL-3R (6H6)-151Eu	Fluidigm Corporation	Cat#3151001B	
Anti-Human CD105/Endoglin (43A3)- 163Dy	Fluidigm Corporation	Cat#3163005B	
Anti-Human CD73 (AD2)-168Er	Fluidigm Corporation	Cat#3168015B	
Anti-Human CD127/IL-7Ra (A019D5)-165Ho	Fluidigm Corporation	Cat#3165008B	
Anti-Human CD197/CCR7 (G043H7)-159Tb	Fluidigm Corporation	Cat#3159003A	
Anti-Human CD163 (GHI/61)-154Sm	Fluidigm Corporation	Cat#3154007B	
Anti-Human CD68 (Y1/82A)-171Yb	Fluidigm Corporation	Cat#3171011B	
Anti-Human CD33 (WM53)-169Tm	Fluidigm Corporation	Cat#3169010B	
Anti-Human CD15/SSEA-1 (W6D3)-164Dy	Fluidigm Corporation	Cat#3164001B	
Anti-Human CD45RA (HI100)-153Eu	Fluidigm Corporation	Cat#3153001B	
Anti-Human CD38 (HIT2)-167Er	Fluidigm Corporation	Cat#3167001B	
	
Biological samples	
	
Human Cells cryopreserved or freshly acquired	Multi-center clinics	N/A	
	
Chemicals, peptides, and recombinant proteins	
	
Maxpar® Cell Acquisition Solution	Fluidigm Corporation	Cat#201241	
Cell-ID™ Cisplatin	Fluidigm Corporation	Cat#201064	
Chromium Next GEM Single Cell 3′ Kit v3.1, 16 rxns	10x Genomics	Cat#1000268	
	
Deposited data	
	
Raw and analyzed data	This paper	GEO: GSE274018
SRA: PRJNA1144164	
Human reference genome NCBI build 38, GRCh38	Genome Reference Consortium	http://www.ncbi.nlm.nih.gov/projects/genome/assembly/grc/human/	
Non-OA Data	Oetjen et al.25	GEO: GSE120221, GSE120446	
	
Experimental models: Organisms/strains	
	
H. sapiens	OA patients from multiple clinical sites	N/A	
	
Software and algorithms	
	
BCL2FASTQ	Illumina	https://support.illumina.com/downloads/bcl2fastq-conversion-software-v2-20.html	
Cellranger V3	10x Genomics	https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/installation	
Seurat V3	Satija Lab	https://satijalab.org/seurat/	
Propeller	GitHub	https://github.com/phipsonlab	
CellChat	cellchatdb	http://www.cellchat.org/	
Custom code	This paper	https://github.com/pchatterjee7/BMAC_OA_nonOA	

Experimental model and study participant details

Ethical approvals

The MILES clinical trial (NCT03818737) was performed in accordance with guidelines and oversight from the Federal Drug Administration (FDA) under IND# 18414 and with subsequent reference to IDE #17894. Study approval was obtained from the Western Institutional Review Board (WIRB) and by Duke and Emory University’s IRB. Only human cells were used for this study.

The following information was supplied relating to ethical approvals: H19004 (GTIRB): Approved 03/12/2019. Based on the Multicenter trial of stem cell therapy for osteoarthritis (MILES). WIRB Protocol #20183019 (Emory University, Duke University, Andrews Institute, Sanford Health (Sioux Falls, Fargo).

Study design and subjects

The non-OA group is of healthy volunteers who were recruited for bone marrow aspiration procedures. The cohort consisted of 10 males and 10 females with ages ranging from 24 to 84 years old and median age of 57 years. A second bone marrow aspiration was performed for 2 donors (Ck, Sk) (biological replicates) either 2 or 5 months after their first aspiration, respectively.25 Mononuclear cells were isolated from bone marrow aspirates via Ficoll density gradient separation, then cryopreserved in 90% FBS/10% DMSO for storage in liquid nitrogen. Each donor’s cryopreserved vials were utilized for assays as detailed in Table S1. For a subset of donors, a second bone marrow aspiration was conducted (for donors Ck and Sk, 2 and 5 months after the initial procedure, respectively) to obtain biological replicates.

The patient data utilized in our study was derived from a large multi-tissue, multi-site clinical trial conducted across the United States. The selection of diverse patients was governed by the trial’s design, which aimed to encapsulate a broad demographic and genetic representation of the population. While the diversity of the patients was determined by the trial’s structure, it nevertheless provided a comprehensive cross-section suitable for observing generalized biological phenomena related to osteoarthritis (OA). Demographic information includes ethnicity and race. Ethnicity is categorized as Hispanic or Latino, Not Hispanic or Latino, Not reported, and Unknown. For racial classification, participants could identify as American Indian, African American/Black, Native Hawaiian/Pacific Islander, Asian, Caucasian/White, or Not reported. These demographic details are critical for interpreting the study results within the context of genetic and cultural backgrounds.

The OA subject group consisted of male and female donors who had been diagnosed with OA and had a Kellgren-Lawrence (KL) scale grade of II, III, or IV, and with no corticosteroid injection for 3 months prior to sample collection and no NSAID use for 6 days prior to sample collection. Donor characteristics and generated data details for both cohorts are summarized in Supplemental Table S1. Samples from the two groups were age and gender matched. To ensure consistency, non-OA samples were cryopreserved using 90% FBS +10% DMSO, while OA samples were processed with Cryostor CS10 containing 10% DMSO.

Method details

Bone marrow aspirate collection and concentration

Bone marrow collected from the posterior superior iliac spine (PSIS) was concentrated by centrifugation (EmCyte GenesisCS Pure BMAC-60 mL) and cryopreserved in Cryostor CS10 (StemCell Technologies) for storage in liquid nitrogen.

Bone marrow aspirate (BMA) from the OA cohort was collected from the posterior superior iliac spine, concentrated by centrifugation (EmCyte GenesisCS Pure BMAC-60 mL) and cryopreserved in Cryostor CS10 (StemCell Technologies) for storage in liquid nitrogen. All BMAC samples were tested for endotoxin and sterility prior to cryopreservation. Assays were performed using matched cryopreserved vials from each donor, collected from the same BMA sample.

Single-cell RNAseq library preparation and sequencing

As an extremely heterogeneous mixture of cells containing diverse immune cell types, the BMAC samples are challenging to process. We performed a pilot study with a few samples to optimize the freezing protocols where we could achieve the highest viability and observed that all the cell types were conserved compared to fresh samples. Upon optimization of the protocol, all samples were cryopreserved in Cryostor10 (Stem cell technologies, CS10) freezing media and were analyzed through droplet-based single-cell RNAseq (10x Genomics Chromium 3′ V3.1). All the cells were thawed for 2 min at 37C and washed and filtered before processing for barcode generation on the 10x Genomics Chromium platform. In case of viability lower than 70%, we used live cell enrichment (Dead cell removal kit, Miltenyi Biotech). The amplified libraries were prepared with ∼5000 cells per sample target. The libraries were sequenced on the Illumina Novaseq6000 platform with an S4 kit to achieve ≥ 30,000 reads per cell.

Single-cell RNA data analysis

The libraries were sequenced on the Illumina sequencing platform, and raw reads were aligned to the human reference genome GRCh38 while cell barcodes were used to assign reads to single cells with CellRanger v3.1. Greater than 230,000 cells were captured and analyzed from both datasets. The diseased datasets were processed with CellRanger software package version 3.1. Our threshold for each sample was set to a minimum cell ≥3, minimum features ≥200, and mitochondrial count ≤40.

After filtering out cells with a low number of genes expressed, high mitochondria gene expressed and low UMIs, 136,484 cells were analyzed in the final analysis. For the non-OA data analysis, we reanalyzed and annotated the published dataset from the raw matrices with the same QC thresholds. The datasets from both diseased and healthy cohorts were initially analyzed with SEURAT version 3, now supplemented with new options available in V4. After filtering, the samples were merged from each group to create OA and Non-OA datasets.

Canonical Correlation Analysis (CCA) was conducted to identify shared sources of variation between conditions/groups, emphasizing the 3000 most variant genes per sample. CCA helps approximate cell alignment based on shared variation. Anchors or mutual nearest neighbors (MNNs) were identified across datasets to provide estimate of batch effect, as described by.73 Incorrectly identified anchors were filtered by assessing the similarity of local neighborhood overlaps. Integration of datasets was achieved using these anchors, transforming cell expression values based on a weighted average determined by cell similarity and anchor scores. Then this object was used as an input for the integration across the conditions. RPCA was used for the large dataset where the datasets are aligned in the PCA space using shared PCs.

Principal components (PCs) of the highly variable genes were computed and the first 30 PCs were included in cluster generation. With a resolution setting of 0.5, 26 cell clusters were detected, and these were visualized in two-dimensional space using uniform manifold approximation and projection (UMAP).74 Cell clusters were distinguished using the Louvain clustering algorithm implemented in Seurat.

After clustering, upon identification of the major cell types to determine the heterogeneity of the BMAC as a whole and, to determine which cell types may differ between the OA and non-OA groups, DEG analysis was performed using the Wilcoxon rank-sum test and MAST to control for batch effects. Once the cell proportions were analyzed, further downstream analysis was performed. The propeller tool75 was used to assess statistical analysis of the cell proportions. Cell communication analysis was performed using CellChat and Nichenet, software applications that predict cellular signals and ligand-target enrichment using network analysis and pattern recognition approaches.26,33 The validation OA dataset of 56 samples was re-analyzed independently using the same QC criteria to have the method reproducibility established. BMI variation analysis was performed on the OA dataset to confirm the effect of BMI on the disease progression.

Flow cytometry

Cryopreserved BMAC samples were prepared and stained with fluorescently conjugated antibodies (Table S4). Samples were first treated with a viability dye (Zombie UV) and incubated with antibody cocktails. Post-staining, samples were fixed using BD Cytofix Fixation Buffer and stored at 4°C. Within three days, samples were processed on a BD LSRFortessa cell analyzer, with over 20,000 events recorded per replicate. Data analysis was performed with FlowJo software, utilizing a standard gating strategy based on scatter, singlets, and viability. The flow cytometry was performed on few of the OA BMAC samples to match the panel of nonOA data. We have later moved toward mass cytometry on all our OA samples to acquire more high throughput data.

Mass cytometry

BMAC samples were readied for mass cytometric staining as per Fluidigm’s instructions. Samples were thawed, washed with serum-free media, and stained with Cell-ID Cisplatin. Staining was neutralized with Maxpar Cell Staining Buffer, followed by cell count and Fc Receptor Blocking Solution incubation. Antibody cocktails were then added. After two washes, samples were fixed with a 1.6% formaldehyde solution, incubated overnight at 4°C in intercalation solution with Cell-ID Intercalator-Ir, washed again, and processed on a Helios mass cytometer using CyTOF software.

Quantification and statistical analysis

Results were analyzed using R v4 and v4.1. All the statistical details of experiments can be found in the STAR Methods section. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001; ns, not significant.

Additional resources

All data are generated from the clinical trial MILES: NCT03818737.

Supplemental information

Document S1. Figures S1‒S12 and Tables S2‒S4

Table S1. Demographic information, related to main Figure 2

A. Donor demographics information Sex, Age, and BMI for OA and nonOA datasets. B. Sample single cell RNAseq generated number of cells, genes from scRNAseq

Supplemental information can be found online at https://doi.org/10.1016/j.isci.2024.110827.
==== Refs
References

1 Glyn-Jones S. Palmer A.J.R. Agricola R. Price A.J. Vincent T.L. Weinans H. Carr A.J. Osteoarthritis Lancet 386 2015 376 387 10.1016/S0140-6736(14)60802-3 25748615
2 Murphy L. Helmick C.G. The impact of osteoarthritis in the United States: a population-health perspective Am. J. Nurs. 112 2012 S13 S19 10.1097/01.NAJ.0000412646.80054.21 22373741
3 Hunter D.J. Bierma-Zeinstra S. Osteoarthritis Lancet 393 2019 1745 1759 10.1016/S0140-6736(19)30417-9 31034380
4 Neogi T. The epidemiology and impact of pain in osteoarthritis Osteoarthritis Cartilage 21 2013 1145 1153 10.1016/j.joca.2013.03.018 23973124
5 Abed E. Bouvard B. Martineau X. Jouzeau J.Y. Reboul P. Lajeunesse D. Elevated hepatocyte growth factor levels in osteoarthritis osteoblasts contribute to their altered response to bone morphogenetic protein-2 and reduced mineralization capacity Bone 75 2015 111 119 10.1016/j.bone.2015.02.001 25667190
6 Blagojevic M. Jinks C. Jeffery A. Jordan K.P. Risk factors for onset of osteoarthritis of the knee in older adults: a systematic review and meta-analysis Osteoarthritis Cartilage 18 2010 24 33 10.1016/j.joca.2009.08.010 19751691
7 Silverwood V. Blagojevic-Bucknall M. Jinks C. Jordan J.L. Protheroe J. Jordan K.P. Current evidence on risk factors for knee osteoarthritis in older adults: a systematic review and meta-analysis Osteoarthritis Cartilage 23 2015 507 515 10.1016/j.joca.2014.11.019 25447976
8 Reijman M. Pols H.A.P. Bergink A.P. Hazes J.M.W. Belo J.N. Lievense A.M. Bierma-Zeinstra S.M.A. Body mass index associated with onset and progression of osteoarthritis of the knee but not of the hip: the Rotterdam Study Ann. Rheum. Dis. 66 2007 158 162 10.1136/ard.2006.053538 16837490
9 Mobasheri A. Henrotin Y. Biomarkers of (osteo)arthritis Biomarkers 20 2015 513 518 10.3109/1354750X.2016.1140930 26954784
10 Bannuru R.R. Osani M.C. Vaysbrot E.E. Arden N.K. Bennell K. Bierma-Zeinstra S.M.A. Kraus V.B. Lohmander L.S. Abbott J.H. Bhandari M. OARSI guidelines for the non-surgical management of knee, hip, and polyarticular osteoarthritis Osteoarthritis Cartilage 27 2019 1578 1589 10.1016/j.joca.2019.06.011 31278997
11 Primorac D. Molnar V. Rod E. Jeleč Ž. Čukelj F. Matišić V. Vrdoljak T. Hudetz D. Hajsok H. Borić I. Knee Osteoarthritis: A Review of Pathogenesis and State-Of-The-Art Non-Operative Therapeutic Considerations Genes 11 2020 854 10.3390/genes11080854 32722615
12 Zhu W. He X. Cheng K. Zhang L. Chen D. Wang X. Qiu G. Cao X. Weng X. Ankylosing spondylitis: etiology, pathogenesis, and treatments Bone Res. 7 2019 22 10.1038/s41413-019-0057-8 31666997
13 Kim G.B. Seo M.S. Park W.T. Lee G.W. Bone Marrow Aspirate Concentrate: Its Uses in Osteoarthritis Int. J. Mol. Sci. 21 2020 3224 10.3390/ijms21093224 32370163
14 Brozovich A. Sinicrope B.J. Bauza G. Niclot F.B. Lintner D. Taraballi F. McCulloch P.C. High Variability of Mesenchymal Stem Cells Obtained via Bone Marrow Aspirate Concentrate Compared With Traditional Bone Marrow Aspiration Technique Orthop. J. Sports Med. 9 2021 23259671211058459 10.1177/23259671211058459
15 Jeyaraman M. Bingi S.K. Muthu S. Jeyaraman N. Packkyarathinam R.P. Ranjan R. Sharma S. Jha S.K. Khanna M. Rajendran S.N.S. Impact of the Process Variables on the Yield of Mesenchymal Stromal Cells from Bone Marrow Aspirate Concentrate Bioengineering (Basel) 9 2022 57 10.3390/bioengineering9020057 35200410
16 Holton J. Imam M. Ward J. Snow M. The Basic Science of Bone Marrow Aspirate Concentrate in Chondral Injuries Orthop. Rev. 8 2016 6659 10.4081/or.2016.6659
17 Imam M.A. Mahmoud S.S.S. Holton J. Abouelmaati D. Elsherbini Y. Snow M. A systematic review of the concept and clinical applications of Bone Marrow Aspirate Concentrate in Orthopaedics SICOT J 3 2017 17 10.1051/sicotj/2017007 29792397
18 Sellam J. Berenbaum F. The role of synovitis in pathophysiology and clinical symptoms of osteoarthritis Nat. Rev. Rheumatol. 6 2010 625 635 10.1038/nrrheum.2010.159 20924410
19 Scanzello C.R. Role of low-grade inflammation in osteoarthritis Curr. Opin. Rheumatol. 29 2017 79 85 10.1097/BOR.0000000000000353 27755180
20 Rosshirt N. Trauth R. Platzer H. Tripel E. Nees T.A. Lorenz H.M. Tretter T. Moradi B. Proinflammatory T cell polarization is already present in patients with early knee osteoarthritis Arthritis Res. Ther. 23 2021 37 10.1186/s13075-020-02410-w 33482899
21 Wang L. He C. Nrf2-mediated anti-inflammatory polarization of macrophages as therapeutic targets for osteoarthritis Front. Immunol. 13 2022 967193 10.3389/fimmu.2022.967193 36032081
22 Kapellos T.S. Bonaguro L. Gemünd I. Human Monocyte Subsets and Phenotypes in Major Chronic Inflammatory Diseases Front. Immunol. 10 2019 2035 10.3389/fimmu.2019.02035
23 Griffin T.M. Scanzello C.R. Innate inflammation and synovial macrophages in osteoarthritis pathophysiology Clin. Exp. Rheumatol. 37 2019 57 63
24 Shi T. Shen X. Gao G. Gene Expression Profiles of Peripheral Blood Monocytes in Osteoarthritis and Analysis of Differentially Expressed Genes BioMed Res. Int. 2019 2019 4291689 10.1155/2019/4291689 31886215
25 Oetjen K.A. Lindblad K.E. Goswami M. Gui G. Dagur P.K. Lai C. Dillon L.W. McCoy J.P. Hourigan C.S. Human bone marrow assessment by single-cell RNA sequencing, mass cytometry, and flow cytometry JCI Insight 3 2018 e124928 10.1172/jci.insight.124928 30518681
26 Jin S. Guerrero-Juarez C.F. Zhang L. Chang I. Ramos R. Kuan C.H. Myung P. Plikus M.V. Nie Q. Inference and analysis of cell-cell communication using CellChat Nat. Commun. 12 2021 1088 10.1038/s41467-021-21246-9 33597522
27 Moazedi-Fuerst F.C. Gruber G. Stradner M.H. Guidolin D. Jones J.C. Bodo K. Wagner K. Peischler D. Krischan V. Weber J. Effect of Laminin-A4 inhibition on cluster formation of human osteoarthritic chondrocytes J. Orthop. Res. 34 2016 419 426 10.1002/jor.23036 26295200
28 Perez-Garcia S. Carrion M. Gutierrez-Canas I. Profile of Matrix-Remodeling Proteinases in Osteoarthritis: Impact of Fibronectin Cells 9 2019 40 10.3390/cells9010040 31877874
29 Lontos K. Adamik J. Tsagianni A. Galson D.L. Chirgwin J.M. Suvannasankha A. The Role of Semaphorin 4D in Bone Remodeling and Cancer Metastasis Front. Endocrinol. 9 2018 322 10.3389/fendo.2018.00322
30 Wang L. Li X. Song Y. Song D. Huang D. The emerging roles of semaphorin4D/CD100 in immunological diseases Biochem. Soc. Trans. 48 2020 2875 2890 10.1042/BST20200821 33258873
31 La Bella S. Rinaldi M. Di Ludovico A. Di Donato G. Di Donato G. Salpietro V. Chiarelli F. Breda L. Genetic Background and Molecular Mechanisms of Juvenile Idiopathic Arthritis Int. J. Mol. Sci. 24 2023 1846 10.3390/ijms24031846 36768167
32 Luo Z. Han Q. Lu J. Ouyang X. Fan Y. Liu Y. Zhou X. Kong J. Liu H. Liu A. Chen D. IL16 Regulates Osteoarthritis Progression as a Target Gene of Novel-miR-81 Cartilage 15 2024 175 183 10.1177/19476035231168387 37086007
33 Browaeys R. Saelens W. Saeys Y. NicheNet: modeling intercellular communication by linking ligands to target genes Nat. Methods 17 2020 159 162 10.1038/s41592-019-0667-5 31819264
34 Upton A.R. Holding C.A. Dharmapatni A.A.S.S.K. Haynes D.R. The expression of RANKL and OPG in the various grades of osteoarthritic cartilage Rheumatol. Int. 32 2012 535 540 10.1007/s00296-010-1733-6 21259010
35 Decout A. Katz J.D. Venkatraman S. Ablasser A. The cGAS-STING pathway as a therapeutic target in inflammatory diseases Nat. Rev. Immunol. 21 2021 548 569 10.1038/s41577-021-00524-z 33833439
36 Thielen N.G.M. van der Kraan P.M. van Caam A.P.M. TGFbeta/BMP Signaling Pathway in Cartilage Homeostasis Cells 8 2019 969 10.3390/cells8090969 31450621
37 Kihara S. Hayashi S. Hashimoto S. Kanzaki N. Takayama K. Matsumoto T. Chinzei N. Iwasa K. Haneda M. Takeuchi K. Cyclin-Dependent Kinase Inhibitor-1-Deficient Mice are Susceptible to Osteoarthritis Associated with Enhanced Inflammation J. Bone Miner. Res. 32 2017 991 1001 10.1002/jbmr.3080 28128866
38 Marinkovic M. Dai Q. Gonzalez A.O. Tran O.N. Block T.J. Harris S.E. Salmon A.B. Yeh C.K. Dean D.D. Chen X.D. Matrix-bound Cyr61/CCN1 is required to retain the properties of the bone marrow mesenchymal stem cell niche but is depleted with aging Matrix Biol. 111 2022 108 132 10.1016/j.matbio.2022.06.004 35752272
39 Zhang Q. Wu J. Cao Q. Xiao L. Wang L. He D. Ouyang G. Lin J. Shen B. Shi Y. A critical role of Cyr61 in interleukin-17-dependent proliferation of fibroblast-like synoviocytes in rheumatoid arthritis Arthritis Rheum. 60 2009 3602 3612 10.1002/art.24999 19950293
40 Mimpen J.Y. Baldwin M.J. Cribbs A.P. Philpott M. Carr A.J. Dakin S.G. Snelling S.J.B. Interleukin-17A Causes Osteoarthritis-Like Transcriptional Changes in Human Osteoarthritis-Derived Chondrocytes and Synovial Fibroblasts In Vitro Front. Immunol. 12 2021 676173 10.3389/fimmu.2021.676173 34054865
41 Barreto G. Manninen M. K Eklund K. Osteoarthritis and Toll-Like Receptors: When Innate Immunity Meets Chondrocyte Apoptosis Biology 9 2020 65 10.3390/biology9040065 32235418
42 Hu X. Li J. Fu M. Zhao X. Wang W. The JAK/STAT signaling pathway: from bench to clinic Signal Transduct. Target. Ther. 6 2021 402 10.1038/s41392-021-00791-1 34824210
43 Ren W. Yin J. Duan J. Liu G. Tan B. Yang G. Wu G. Bazer F.W. Peng Y. Yin Y. mTORC1 signaling and IL-17 expression: Defining pathways and possible therapeutic targets Eur. J. Immunol. 46 2016 291 299 10.1002/eji.201545886 26558536
44 Guo Q. Chen X. Chen J. Zheng G. Xie C. Wu H. Miao Z. Lin Y. Wang X. Gao W. STING promotes senescence, apoptosis, and extracellular matrix degradation in osteoarthritis via the NF-kappaB signaling pathway Cell Death Dis. 12 2021 13 10.1038/s41419-020-03341-9 33414452
45 Hari P. Millar F.R. Tarrats N. Birch J. Quintanilla A. Rink C.J. Fernández-Duran I. Muir M. Finch A.J. Brunton V.G. The innate immune sensor Toll-like receptor 2 controls the senescence-associated secretory phenotype Sci. Adv. 5 2019 eaaw0254 10.1126/sciadv.aaw0254 31183403
46 Schubert M. Klinger B. Klünemann M. Sieber A. Uhlitz F. Sauer S. Garnett M.J. Blüthgen N. Saez-Rodriguez J. Perturbation-response genes reveal signaling footprints in cancer gene expression Nat. Commun. 9 2018 20 10.1038/s41467-017-02391-6 29295995
47 Katsoula G. Kreitmaier P. Zeggini E. Insights into the molecular landscape of osteoarthritis in human tissues Curr. Opin. Rheumatol. 34 2022 79 90 10.1097/BOR.0000000000000853 34750308
48 Haringman J.J. Smeets T.J.M. Reinders-Blankert P. Tak P.P. Chemokine and chemokine receptor expression in paired peripheral blood mononuclear cells and synovial tissue of patients with rheumatoid arthritis, osteoarthritis, and reactive arthritis Ann. Rheum. Dis. 65 2006 294 300 10.1136/ard.2005.037176 16107514
49 Zhao E. Xu H. Wang L. Kryczek I. Wu K. Hu Y. Wang G. Zou W. Bone marrow and the control of immunity Cell. Mol. Immunol. 9 2012 11 19 10.1038/cmi.2011.47 22020068
50 Covre L.P. De Maeyer R.P.H. Gomes D.C.O. Akbar A.N. The role of senescent T cells in immunopathology Aging Cell 19 2020 e13272 10.1111/acel.13272 33166035
51 Faust H.J. Zhang H. Han J. Wolf M.T. Jeon O.H. Sadtler K. Peña A.N. Chung L. Maestas D.R. Jr. Tam A.J. IL-17 and immunologically induced senescence regulate response to injury in osteoarthritis J. Clin. Invest. 130 2020 5493 5507 10.1172/JCI134091 32955487
52 Boer C.G. Hatzikotoulas K. Southam L. Stefánsdóttir L. Zhang Y. Coutinho de Almeida R. Wu T.T. Zheng J. Hartley A. Teder-Laving M. Deciphering osteoarthritis genetics across 826,690 individuals from 9 populations Cell 184 2021 4784 4818.e17 10.1016/j.cell.2021.07.038 34450027
53 Sun K. Chen X.S. Muzhylko T. Andrade F.C.D. Doctors' recommendations and healthy lifestyle behaviors among individuals with hypertension in Brazil Prev. Med. Rep. 35 2023 102315 10.1016/j.pmedr.2023.102315 37576845
54 Li Y.S. Luo W. Zhu S.A. Lei G.H. T Cells in Osteoarthritis: Alterations and Beyond Front. Immunol. 8 2017 356 10.3389/fimmu.2017.00356 28424692
55 Keindl M. Davies R. Bergum B. Brun J.G. Hammenfors D. Jonsson R. Lyssenko V. Appel S. Impaired activation of STAT5 upon IL-2 stimulation in Tregs and elevated sIL-2R in Sjogren's syndrome Arthritis Res. Ther. 24 2022 101 10.1186/s13075-022-02769-y 35526080
56 Prockop D.J. Oh J.Y. Mesenchymal stem/stromal cells (MSCs): role as guardians of inflammation Mol. Ther. 20 2012 14 20 10.1038/mt.2011.211 22008910
57 Zou M.L. Chen Z.H. Teng Y.Y. Liu S.Y. Jia Y. Zhang K.W. Sun Z.L. Wu J.J. Yuan Z.D. Feng Y. The Smad Dependent TGF-beta and BMP Signaling Pathway in Bone Remodeling and Therapies Front. Mol. Biosci. 8 2021 593310 10.3389/fmolb.2021.593310 34026818
58 Lee H.S. Millward-Sadler S.J. Wright M.O. Nuki G. Al-Jamal R. Salter D.M. Activation of Integrin-RACK1/PKCalpha signalling in human articular chondrocyte mechanotransduction Osteoarthritis Cartilage 10 2002 890 897 10.1053/joca.2002.0842 12435334
59 Xie W. Xiao W. Tang K. Zhang L. Li Y. Yes-Associated Protein 1: Role and Treatment Prospects in Orthopedic Degenerative Diseases Front. Cell Dev. Biol. 8 2020 573455 10.3389/fcell.2020.573455 33178690
60 Liu T. Zhang L. Joo D. Sun S.C. NF-kappaB signaling in inflammation Signal Transduct. Target. Ther. 2 2017 17023 10.1038/sigtrans.2017.23 29158945
61 Liu Z. Wang P. Cen S. Gao L. Xie Z. Wu X. Su H. Wu Y. Shen H. Increased BMPR1A Expression Enhances the Adipogenic Differentiation of Mesenchymal Stem Cells in Patients with Ankylosing Spondylitis Stem Cells Int. 2019 2019 4143167 10.1155/2019/4143167 31827527
62 Benigni G. Dimitrova P. Antonangeli F. Sanseviero E. Milanova V. Blom A. van Lent P. Morrone S. Santoni A. Bernardini G. CXCR3/CXCL10 Axis Regulates Neutrophil-NK Cell Cross-Talk Determining the Severity of Experimental Osteoarthritis J. Immunol. 198 2017 2115 2124 10.4049/jimmunol.1601359 28108560
63 Christopoulos P.F. Gjolberg T.T. Kruger S. Haraldsen G. Andersen J.T. Sundlisaeter E. Targeting the Notch Signaling Pathway in Chronic Inflammatory Diseases Front. Immunol. 12 2021 668207 10.3389/fimmu.2021.668207 33912195
64 Hou S.M. Chen P.C. Lin C.M. Fang M.L. Chi M.C. Liu J.F. CXCL1 contributes to IL-6 expression in osteoarthritis and rheumatoid arthritis synovial fibroblasts by CXCR2, c-Raf, MAPK, and AP-1 pathway Arthritis Res. Ther. 22 2020 251 10.1186/s13075-020-02331-8 33087182
65 Huang Y. Huang Q. Su H. Mai X. Feng E. Cao Z. Zeng X. TAR DNA-binding protein 43 inhibits inflammatory response and protects chondrocyte function by modulating RACK1 expression in osteoarthritis Biomed. Pharmacother. 85 2017 362 371 10.1016/j.biopha.2016.11.037 27890430
66 Kam N.W. Liu D. Cai Z. Mak W.Y. Wong C.K. Chiu K.H. Wong K.Y. Tsang W.L. Tam L.S. Synoviocytes-derived Interleukin 35 Potentiates B Cell Response in Patients with Osteoarthritis and Rheumatoid Arthritis J. Rheumatol. 45 2018 563 573 10.3899/jrheum.161363 29247146
67 Martinez-Anton A. Gras D. Bourdin A. Dubreuil P. Chanez P. KIT as a therapeutic target for non-oncological diseases Pharmacol. Ther. 197 2019 11 37 10.1016/j.pharmthera.2018.12.008 30557630
68 Merino A. Buendia P. Martin-Malo A. Aljama P. Ramirez R. Carracedo J. Senescent CD14+CD16+ monocytes exhibit proinflammatory and proatherosclerotic activity J. Immunol. 186 2011 1809 1815 10.4049/jimmunol.1001866 21191073
69 Hombach A.A. Geumann U. Gunther C. Hermann F.G. Abken H. IL7-IL12 Engineered Mesenchymal Stem Cells (MSCs) Improve A CAR T Cell Attack Against Colorectal Cancer Cells Cells 9 2020 873 10.3390/cells9040873 32260097
70 Han Y. Yang J. Fang J. Zhou Y. Candi E. Wang J. Hua D. Shao C. Shi Y. The secretion profile of mesenchymal stem cells and potential applications in treating human diseases Signal Transduct. Target. Ther. 7 2022 92 10.1038/s41392-022-00932-0 35314676
71 Khalid R.S. Khan I. Zaidi M.B. Naeem N. Haneef K. Qazi R.E.M. Habib R. Malick T.S. Ali A. Salim A. IL-7 overexpression enhances therapeutic potential of rat bone marrow mesenchymal stem cells for diabetic wounds Wound Repair Regen. 27 2019 235 248 10.1111/wrr.12706 30761686
72 Shao Y. Zhao C. Pan J. Zeng C. Zhang H. Liu L. Fan K. Liu X. Luo B. Fang H. BMP5 silencing inhibits chondrocyte senescence and apoptosis as well as osteoarthritis progression in mice Aging (Albany NY) 13 2021 9646 9664 10.18632/aging.202708 33744859
73 Stuart T. Butler A. Hoffman P. Hafemeister C. Papalexi E. Mauck W.M. 3rd Hao Y. Stoeckius M. Smibert P. Satija R. Comprehensive Integration of Single-Cell Data Cell 177 2019 1888 1902.e21 10.1016/j.cell.2019.05.031 31178118
74 Becht E. McInnes L. Healy J. Dutertre C.A. Kwok I.W.H. Ng L.G. Ginhoux F. Newell E.W. Dimensionality reduction for visualizing single-cell data using UMAP Nat. Biotechnol. 37 2018 38 44 10.1038/nbt.4314
75 Phipson B. Sim C.B. Porrello E.R. Hewitt A.W. Powell J. Oshlack A. propeller: testing for differences in cell type proportions in single cell data Bioinformatics 38 2022 4720 4726 10.1093/bioinformatics/btac582 36005887
