Single-Cell Transcriptomics of Multi-Site Cell Therapy in Osteoarthritis: Tissue-Specific Traits and Treatment Correlations

Single-Cell Transcriptomics of Multi-Site Cell Therapy in Osteoarthritis: Tissue-Specific Traits and Treatment Correlations | medRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-P4HH5NV'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search Single-Cell Transcriptomics of Multi-Site Cell Therapy in Osteoarthritis: Tissue-Specific Traits and Treatment Correlations View ORCID Profile Paramita Chatterjee , Hazel Y. Stevens , Linda E. Kippner , Carolyn Yeago , Hicham Drissi , Kenneth Mautner , Scott Boden , View ORCID Profile Greg Gibson , View ORCID Profile Krishnendu Roy doi: https://doi.org/10.1101/2025.04.14.25325743 Paramita Chatterjee 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 Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Paramita Chatterjee Hazel Y. Stevens 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 Find this author on Google Scholar Find this author on PubMed Search for this author on this site Linda E. Kippner 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 Find this author on Google Scholar Find this author on PubMed Search for this author on this site Carolyn Yeago 2 The Parker H. Petit Institute for Bioengineering and Biosciences Georgia Institute of Technology , Atlanta, GA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Hicham Drissi 3 Department of Orthopedics, Emory University School of Medicine , Atlanta, GA, 30322, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Kenneth Mautner 3 Department of Orthopedics, Emory University School of Medicine , Atlanta, GA, 30322, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Scott Boden 3 Department of Orthopedics, Emory University School of Medicine , Atlanta, GA, 30322, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Greg Gibson 4 School of Biological Sciences, Georgia Institute of Technology , Atlanta GA 30332, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Greg Gibson Krishnendu Roy 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 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 Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Krishnendu Roy For correspondence: krish.roy{at}vanderbilt.edu Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract Knee-Osteoarthritis (Knee-OA) is a prevalent joint disorder lacking FDA-approved cell-therapies to halt its progression. This study uses single-cell RNA sequencing to analyze bone marrow aspirate concentrate (BMAC) and stromal vascular fraction (SVF) samples in a clinical trial of autologous cell therapies. Trial site-specific variability was significant in BMAC, necessitating tailored normalization, whereas SVF was less affected, likely due to uniform subcutaneous fat sampling. Variance partitioning and tensor decomposition identified site effects in BMAC but revealed shared pathways across cell types in both tissues. Differential gene expression (DEG) analysis between responders and non-responders yielded no significant findings, though likelihood ratio testing (LRT) revealed enrichment for DEG patterns linked to disease severity, potentially masked by patient heterogeneity. Key BMAC pathways included oxidative phosphorylation, unfolded protein response, and TNFα signaling. Cell-cell communication analysis suggested enhanced HLA signaling in non-responder MSCs, consistent with inflammation, while responders showed more coordinated immune interactions. BMAC-MSCs promoted chondrocyte proliferation, whereas SVF-MSCs emphasized immune regulation. Variability in therapy outcomes reflects patient heterogeneity beyond genomic factors, complicating the immediate use of genomic profiling to guide treatment. Nonetheless, as molecular pathways are better understood, integrating genomic insights into personalized strategies may become feasible. Introduction Osteoarthritis (OA) is a degenerative joint disease that significantly diminishes the quality of life for millions of people worldwide. It is characterized primarily by debilitating joint pain and progressive cartilage loss 1 . Traditional treatment modalities, such as physical therapy, pharmacological interventions, and intra-articular corticosteroid injections, provide symptomatic relief but often fail to arrest disease progression or restore damaged tissues effectively. In recent years, the emergence of regenerative medicine has improved the prospective therapeutic landscape for OA by introducing innovative cellular therapies designed to exploit the body’s inherent regenerative capabilities 2 . Among these therapies, autologous bone marrow aspirate concentrate (BMAC) 3 and adipose-derived stromal vascular fraction (SVF) 4 have gained prominence. BMAC is rich in stromal cells that aim to harness the regenerative and immunomodulatory properties of bone marrow-derived cells to tackle the multifaceted pathology of OA, which includes cartilage degradation, subchondral bone remodeling, and inflammation. Similarly, SVF 5 , derived from adipose tissue, contains a heterogeneous mix of cells with potent anti-inflammatory and regenerative properties, offering a promising avenue for managing knee OA by mitigating inflammation and promoting repair within the joint. Despite the promise and growing clinical application of these therapies, a thorough evaluation of their efficacy, particularly in comparison to other cellular therapies and standard of care, remains a critical unmet need 6 . To address this, a comprehensive phase 2/3 multicenter, single-blind, randomized controlled clinical trial, MILES ( NCT03818737 ), was undertaken. This trial aimed to assess the safety and efficacy of BMAC, SVF, and allogeneic human umbilical cord tissue-derived mesenchymal stromal cells against conventional treatments for knee OA in a robust study design involving 480 patients 7 . The trial’s methodology involved randomized allocation of patients into treatment arms receiving cellular therapies or corticosteroid treatments, with follow-up extending to three years. Central to this study was the innovative integration of tissue-specific single-cell RNA sequencing analysis, which was employed to probe the molecular and cellular dynamics within bone marrow and adipose. Transcriptomic profiling was performed to determine whether there is a correlation between the cellular dynamics observed in treatments like BMAC and clinical outcomes such as pain relief and improved knee function. Association between molecular profiles and patient-reported outcomes enhances our understanding of how cellular therapies modulate the OA environment at a molecular level, potentially leading to more tailored and effective treatment strategies. Our previous study 8 reported on the use of single-cell transcriptomics to analyze bone marrow aspirate concentrate (BMAC) from OA patients, revealing significant alterations in key cell types, including mesenchymal stromal cells (MSCs), natural killer (NK) cells, monocytes, B cells, and T cells between OA patients and controls. The findings suggested that bone marrow plays a significant role as a reservoir for cellular mediators of inflammation and immune modulation in OA. The current study expands on this work by incorporating additional BMAC profiles along with data from stromal vascular fractions (SVF) from OA patients and exploring their relationship to therapeutic response. This study thus provides a more detailed characterization of the cellular environments and interactions within these tissues. Overall, the inter-individual variability of profiles obscures any major signatures of remission or disease progression, though we do identify gene expression features that may be enriched in some recipients with poor outcomes. Efforts to incorporate personalized genomic profiling into therapeutic decision making will need to account for the patient variability. Materials and Methods Study Design and Participants The MILES 9 study was a multicenter single-blinded clinical trial performed in accordance with guidelines and oversight from the FDA, and under the management of a contracted research organization. Study approval was obtained from the Western Institutional Review Board and by Duke and Emory University’s institutional review boards. In accordance with the FDA requirements, an investigational new drug application was filed (#18414), referencing investigative device exemption #17894. The study was also registered with ClinicalTrials.gov 7 ( NCT03818737 ). A total of 570 patients were screened to identify 480 eligible patients who were randomized at five clinical sites across five states within the United States. The five clinical sites included clinics from Emory University (GA), Sanford Health (SD and ND), the Andrews Institute (FL), and Duke University (NC). Participants were randomized to a four-cohort parallel-design study. Participants were further subdivided into three arms, each of which was matched with one third (n=40) of the corticosteroid injection (CSI) controls ( Figure 1A ), arm 1—autologous bone marrow concentrate (BMAC) (n = 120); arm 2—autologous SVF (n = 120); arm 3—allogeneic umbilical cord tissue (UCT) (n = 120). The control cohort of CSI was subsequently aggregated to allow a 1:1:1:1 comparison for analysis. Subjects were enrolled if they were between 40 and 70 years of age and carried a diagnosis of knee OA as determined by radiographs within 3 months of their clinical visit, also having a Kellgren-Lawrence (KL) OA severity score between 2-4. A full list of eligibility criteria is provided on ClinicalTrials.gov ( NCT03818737 ). In the BMAC arm, a total of 12 subjects did not complete the study: 1 related to release criteria, 10 related to subject withdrawal, and 1 lost to follow-up. In the SVF cohort, 12 subjects did not complete the study: 1 due to investigator withdrawal, 4 related to release criteria, 6 related to subject withdrawal, and 1 due to screen failure. In the UCT cohort, 13 subjects did not complete the study: 2 due to investigator withdrawal, 10 related to subject withdrawal, and 1 lost to follow-up. In the CSI cohort, 11 subjects did not complete the study: 1 due to investigator decision, 6 related to subject withdrawal, and 4 lost to follow-up. Of the remaining patients, samples from 134 patients from each of the BMAC and SVF arms, and samples from all 8 UCT-MSC lots met minimum criteria for additional processing and were analyzed for single cell RNAseq ( Figure 1B ). Patients returned to the clinic for MRI to assess cartilage and joint health 10 at baseline, 6 months, and 1 year. Download figure Open in new tab Figure 1: Experimental and Analytical Design and for scRNAseq data A. MILES Clinical Trial Data Collection Design: The trial consisted of 480 total participants, randomized into one of 3 arms. In each arm, patients received either cell therapy injection or corticosteroid injection (control) to the knee. BMAC and SVF cell product consisted of autologous therapeutic cells collected from the patient’s bone marrow or adipose tissue, respectively. The UCT-MSC allogeneic cell therapy was manufactured from umbilical cord tissue, with 8 lots included. B . Out of the total enrolled patients in the study, samples from 134 patients from each of the BMAC and SVF arms, and 8 samples from UCT-MSC lots were analyzed using 10x Genomics 3’ scRNAseq kit V3.1. All the autologous data from BMAC and SVF are collected pre-treatment. For BMAC, among 134, there are 32 patients treated with corticosteroid injection (CSI), 11 patients dropped out, 29 patients were grouped as non-responders as of VAS responsiveness and 62 patients were responders as of the VAS responsiveness after 12 months follow-up. For SVF, among 134, there are 24 patients treated with corticosteroid injection (CSI), 30 patients dropped out, 33 patients were grouped as non-responders as of VAS responsiveness and 47 patients were responders as of the VAS responsiveness after 12 months follow-up. For UCT-MSCs being allogeneic, there were data generated from 8 donor samples who are not OA patients. C. Data analysis pipeline design Randomization and Masking Subjects were stratified by clinical center, and after informed consent was obtained and eligibility criteria verified, they were randomized to treatment groups 11 . Treatment assignments were also stratified by clinical center and generated using a pseudo-random number generator with randomly permuted blocks. These assignments were stored in the Medidata Rave cloud-based data management system developed by the contracted research organization. All subjects underwent the harvesting procedure per their assigned study arm, and then, according to the randomization scheme, they were blinded to the actual treatment received (e.g., SVF versus CSI, BMAC versus CSI, or cells versus CSI). As a single-blinded study, the site principal investigators were not required to be blinded; however, physicians and subjects were blinded to the injection. The blinding was implemented by limiting visualization of the syringe contents with opaque coverings, in addition to performing sham BMAC and SVF harvests in patients within certain cohorts. Single-Cell RNA Sequencing Analysis To promote consistency, standardized protocols were followed across the five study sites for cellular harvesting, product preparation, injection procedures, and corticosteroid injections (CSI). While BMAC and SVF samples were fresh autologous products, MSCs derived from umbilical cord tissue had been cryopreserved and prepared to cGMP standards at Duke University. In our earlier publication 8 , 77 BMAC samples from the MILES cohort were used to compare healthy and disease conditions; an additional 57 samples in this study expanded the total to 134, following the same protocols. Single-cell suspensions from cryopreserved BMAC and SVF samples were prepared and assessed using the NC-200 Nucleocounter cell counter. Transcriptome analysis was conducted on both the BMAC and SVF cells, as well as the allogeneic UCT-MSC samples. Library preparation and sequencing were performed using the 10x Genomics Chromium platform with the 3′ V3.1 kit on an Illumina Novaseq 6000. The analysis pipeline is depicted in Figure 1C . Raw data were processed with Cell Ranger (10X Genomics) to generate count matrices, which were then analyzed with the Seurat (v4.0) R package 12 . Cells were filtered by quality, requiring 500-10,000 UMIs per cell and <20% mitochondrial gene content, with genes detected in fewer than 3 cells excluded. For BMAC samples, the median cell count per sample was 3680 (s.d. 2218) and the median reads per cell were 45282 (s.d. 86678). For SVF samples, the median cell count per sample was 4869 (s.d. 2058) and the median reads per cell were 105555 (s.d. 3746). After normalization, principal component analysis (PCA) was used to reduce dimensionality. Among several normalization approaches tested, including robust principal component analysis (rPCA) 13 , canonical correlation analysis (CCA) and Harmony 14 , rPCA was selected for its ability to balance batch effects with biological signal retention. Clustering was performed using the Louvain algorithm, and after quality control comparison we settled on parameters that identified 37 cell clusters. Cell Type Annotation and Abundance Assessment Identified clusters were annotated with cell types based on the expression of known marker genes from literature and databases. The DISCO 15 database ( https://www.immunesinglecell.org/ct_pub ) was utilized for automated annotation, supplemented by manual review. For ambiguous clusters, additional marker gene analyses and literature review were performed to assign the most probable cell type. The results were visualized using UMAP 16 , highlighting the distribution and relationship of the annotated cell types within the samples. A matrix of cell counts per identified cell type for each sample was prepared by aggregating the UMI counts for all genes within a cell, thereby assigning each cell to its annotated cell type based on the clustering and annotation steps previously described. Subsequently, the relative abundance of each cell type across different samples and conditions was quantified using the Propeller tool 17 . Pseudobulk Differential Gene Expression Analysis To investigate the relationship between differential gene expression and clinical response in the trial arms, pseudobulk 18 profiles were generated with customized scripts by aggregating single-cell RNA-seq counts across annotated cell types within each sample. Initial quality control included principal component analysis (PCA) and hierarchical clustering on log-transformed or variance-stabilized data to assess sample similarity and detect outliers. To mitigate site-specific effects across clinical centers, batch correction was applied to the pseudobulk data using CombatSeq 19 . Identification of differentially expressed genes (DEGs) was then performed using DESeq2 20 . For BMAC, we utilized likelihood ratio tests (LRT) 21 to compare full and reduced models (excluding age, sex, and BMI) to identify enriched sets of DEGs. In parallel, variance component partitioning was performed with dreamlet 22 , 23 . For the SVF dataset, in addition to the pseudobulk analyses, single-cell RNA-seq data was also analyzed with scITD which uses Tucker tensor decomposition to capture gene sets co-regulated across multiple cell types. All results were adjusted for multiple comparisons using the Benjamini-Hochberg method with a 5% false discovery rate (FDR) 24 threshold. Cell-to-Cell Communication Analysis CellChat 25 and NicheNet 26 were used to identify likely signaling networks and pathways that may mediate intercellular interactions. CellChat was used to analyze and visualize cell-cell communication networks based on inferred ligand-receptor interactions previously annotated in a comprehensive well-validated database. NicheNet further generates sender-receiver models of cellular communication, prioritizing ligand-target links that are most likely to explain observed gene expression changes in receiver cells. Both tools were implemented with standard workflows recommended by the developers. Pathway Enrichment Analysis Pathway enrichment analysis was conducted using a combination of Gene Set Enrichment Analysis (GSEA) 27 and its optimized counterpart, fast Gene Set Enrichment Analysis (fgsea) 28 , to rank genes based on their differential expression and assess the enrichment of predefined gene sets. The ESCAPE R 29 package was employed for the integration of pathway databases, enabling comprehensive visualization and interpretation of enriched pathways. Statistical significance was determined using nominal p-values and adjusted for multiple comparisons also with a 5% FDR threshold. Results Cluster Analysis Reveals Distinct Cell Type Compositions in BMAC and SVF Datasets Cluster analysis of the BMAC and SVF datasets revealed 18 and 19 distinct cell-type clusters, respectively ( Figure 2A and 2B ). In the BMAC dataset, major cell types identified included immune cells, monocytes, erythrocytes, megakaryocytes (MEGA), mesenchymal stromal cells (MSCs), and hematopoietic stem cells (HSPCs), with MSCs comprising a low percentage, approximately 0.2% of the total cell population after filtering and quality control. In the SVF dataset, predominant cell types included adipose-derived stromal cells, epithelial cells, fibroblasts, macrophages, and MSC which constituted approximately 5% of the total cell population. The distinction between adipose-derived stromal cells (ASCs) and MSCs was based on marker expression: ASCs exhibited high CD34 and HLA-DR expression but low CD90, while MSCs were characterized by high CD90 and low HLA-DR expression. Following extensive comparison of different procedures for integration of datasets obtained from the different collection sites, we settled on rPCA for all subsequent analyses. As with all single cell RNAseq analyses, the final cluster assignment represents a balance between consistency across datasets, noise reduction, and matching cell types to biological expectation. Erythrocytes were removed from the final clustering as their proportion is highly variable and obscures trends in proportions of the other cell types. Download figure Open in new tab Figure 2: UMAPs and Cellular compositions for scRNAseq celltype clusters for BMAC and SVF datasets A . Celltypes identified for BMAC datasets showing cellular composition. 18 large celltype clusters are identified. The erythrocyte cluster was excluded to reduce the noise for BMAC. Large immune celltypes and monocytes are identified. B . BMAC cell composition barplot aligned with clinical collection site. B . Celltype identified for SVF datasets showing the cellular composition of SVF as a product. 19 large celltype clusters are identified among these immune cell, epithelial cell, macrophages, adipose derived stromal cells and fibroblasts are identified. D . SVF cell composition barplot aligned with clinical collection site. Dittoseq package is used for the barplots, and Propeller tool is used for the statistical analysis on the cell proportions to investigate any significant correlation between cell types and the metadata like Sex, Site, Age, Individuals, and Treatment responses. Cell Abundance Analysis Inferred cell type proportions per donor group by collection site are shown in Figure 2C, D for BMAC and SVF respectively. Using the Propellor tool to evaluate variability of cell type abudance, in the BMAC, macrophages, progenitor B and plasma cells, CD14 monocytes, and megakaryocytes were all found to vary among sites regardless of normalization method ( Table S1 ). This is most likely due to subtle differences in bone marrow sampling by different physicians. Since pro-B and mature B cells, dendritic cells, MSC as well as naïve CD4 and CD8 T cells were additionally significantly variable among sites when using Harmony, we regarded this method as less robust but show the proportions in Figure S1 . In contrast, variability across normalization methods was not observed for the SVF datasets shown in Figure 2D . After controlling for collection site, age, sex, BMI and Kellgren-Lawrence disease severity index, no cell type proportions varied significantly between responders and non-responders ( Tables S1-S4 ). Variance Partition Analysis for BMAC and SVF Turning our attention to gene expression within cell types, dreamlet was used to perform variance partition analysis of the influence of collection site ( Figure 3A ) and therapeutic response status ( Figure 3B ). The violin plots show the range of variance explained by these factors for each expressed gene with the main effect shown in brown and residual unexplained variance in gray. These results imply a substantial influence of collection site on BMAC gene expression variability, with a median between 10% and 25% in most BMAC cell types, for some genes over 50% of the variability due to site, most notably in MSCs and plasma B cells. Only CD16 positive monocytes showed an elevated collection site influence in the SVF dataset. Download figure Open in new tab Figure 3: Assessment of site variance impact on gene expression in two separate datasets. A . This figure depicts the impact of the site variable on gene expression for various cell types within two independent datasets, visualized side by side for comparison (Left: BMAC, Right: SVF). Violin plots are used to represent the distribution of variance percentages explained by the site variable within each cell type. This plot illustrates BMAC with cell types from CD14 Monocytes to Transitional B cells. The black dots indicate the median variance explained by the site, and the vertical lines delineate the interquartile range for each cell type. The distribution and spread of the violin plots convey the degree to which the site variable influences gene expression within each cell type. A broader plot suggests a greater diversity in the variance explained by the site, whereas a narrower plot indicates less variability. By comparing the range and median values of the variance explained across the two panels, we can infer which dataset exhibits a more substantial effect of the site on gene expression. This analysis employs the VariancePartition from Dreamlet package, which utilizes linear mixed models to quantify the fractional contribution of the site variable to the overall expression variation. In strong contrast, treatment response contributed minimally to the variance in both BMAC and SVF datasets, explaining less than 10-15% of the variability of all genes across most cell types. Interestingly, the only exception was MSCs in the BMAC with approximately 20% of the variance in gene expression associated with treatment response, though this result is tempered by the relatively small counts of MSC in that tissue compartment. There was also a slightly elevated response of ASC/MSC and CD16 monocyte gene expression in the SVF. Identification of Multicellular Gene Expression Programs Across Diverse Cell Types Using scITD in Osteoarthritis SVF: Insights into Donor-Specific Variabilities and Immune-Stromal Interactions We next performed tensor decomposition to explore how gene expression programs shared across multiple cell types associate with biological or technical aspects of the study. For this, we used scITD software, which provided the most informative results for gene expression programs associated with immune and stromal cell interactions in the SVF dataset. The software first identifies components of variation shared across cell types and shows how these covary with collection site or treatment response ( Figure 4A ) and then allows visualization of the expression of key transcripts associated with the components in the contributing cell types. Factors 2 and 4 weakly associated with site, but none of the first five factors associated with response. Nevertheless, Factor 1 exhibited a strong association with specific immune cell types ( Figure 4B ), notably antibody-secreting cells, B cells, CD4 T cells, NK cells, macrophages, and monocytes. Upregulated genes associated with Factor 1 included FBN1, IGFBP5 , and MFAP5 , all potentially engaged in immune-stromal interactions, while TGFBR3 and PRG4 were down-regulated in these cell types. Gene Set Enrichment Analysis (GSEA) of the Factor 1 loadings further supported stromal cell enrichment for pathways related to innate immune response, B cell-mediated immunity, and T cell activation. In parallel, pathways linked to cytokine signaling, antigen presentation, and immune effector functions were highlighted in the immune cells in the SVF. Loadings of other gene sets and cell types for Factors 2 through 5 are shown in Figure S2 . Download figure Open in new tab Figure 4: Tensor Factorization of Gene Expression Profiling in SVF. A . The heatmap depicts z-score normalized expression data across multiple donors, stratified by site and condition, with sidebars indicating sample origin (Site) and clinical outcome (Condition). Red shades indicate higher factor levels, and blue shades represent lower factor levels relative to the mean. B . The heatmap illustrates the association of various genes with ‘Factor 1’ across different cell types. Red intensity signifies a stronger association. Listed genes demonstrate a pronounced association with Factor 1, suggesting a potential biological significance in these specific cell types Pseudobulk Analysis of BMAC shows enrichment of DEGs correlated with Clinical Response Reflecting the above documented strong site effect, scITD Factor 1 for BMAC was strongly associated with collection site. Since the remaining BMAC factors did not clearly associate with a biological covariate, an alternative analytical approach was taken to better address whether response to therapy correlates with differential gene expression. Single cell profiles were converted to pseudobulk profiles for each cell type in each donor. Since there was no significant differential gene expression between responders and non-responders below the 5% false discovery rate threshold in any of the individual cell types after standard hypothesis testing, we utilized the likelihood ratio test (LRT) option in DESeq2 to contrast a null model of no differences among response groups with alternate models. We thus searched for genes whose profile in pseudobulk cells representing each cell type in all individuals was more likely to fit a specific pattern of differential expression with respect to treatment response. Figure 5A shows that in all 18 BMAC cell types by far the most commonly supported model showed a gradient from lowest expression in samples from individuals who did not complete the study to highest in responders, who generally resemble control. The number of genes following this pattern ranged from 2,923 in mature B cells (and similarly in the major T cell sub-types) to fewer than 100 in MSC and macrophages (reflecting the small number of such cells). Download figure Open in new tab Figure 5: Differential gene expression analysis across treatment groups using DESeq2 likelihood ratio test (LRT). A. The plots visualize the z-score of gene abundance for multiple cell type containing groups of DEGs with a consistent pattern. Each subplot represents a distinct group, labeled by group number and the count of genes analyzed. The categories on the x-axis denote the different treatment conditions: control, incomplete responder, non-responder, and responder. Points represent individual genes, and the lines connect the median z-scores across the treatment conditions, showing the relative expression levels. Box plots overlay the data points to summarize the distribution of expression within each condition. The color red signifies genes belonging to ‘one’ category, as indicated by the color legend. B . The pattern is consistent after removing the site variability from the data. C . Key pathways identified in heatmap showcases various inflammatory and degenerative differential pathways, with significance denoted by asterisks in the non_responders. GSEA database is used for the pathway analysis. Further adjustment was performed using Combat to remove the batch effect of collection site on the pseudobulk expression profiles. This considerably reduced the numbers of significant genes as much as two orders of magnitude, yet around 100 genes with this pattern were retained, with examples show in Figure 5B . Gene set enrichment analysis (GSEA) of these DEGs identified key biological pathways impacted by the treatment. Upregulated pathways included KRAS signaling and hypoxia, while downregulated pathways encompassed oxidative phosphorylation, unfolded protein response, and TNFα signaling via NFκB ( Figure 5C ). Restricting the analysis to the responder and non-responder categories shows that the lack of significant differential gene expression between groups may be due to the high variability within groups. This is clearly illustrated by the dot plots for the ten most differentially expressed genes in four cell types in Figure 6 (need to decide if we want to move to supplement or keep in main figures) . The size of each dot represents the proportion of cells for the donor who express the transcript, and the shading represents the average gene expression. We chose ten individuals of each class at random from across the collection sites, and grouped them by response group, but visually it is clear that there is little correspondence between response status and expression. The four selected cell types—CD16 monocytes, Pro-B cells, Plasma B cells, and MSCs—display variation in gene expression without clear separation between responders and non-responders. In CD16 monocytes, for example, while most of the ten genes are co-expressed at a higher level in non-responders, one responder (labelled 04bmac039) is among the high-expressors, while others have elevated expression of one or a few of the genes. A similar pattern by the other cell types, with the MSCs showing highly variable expression of multiple ribosomal protein-encoding genes. These results indicate that gene expression variation does not distinctly separate responders from non-responders. Download figure Open in new tab Figure 6: Expression Variability of Selected Genes in BMAC Cell Types Across Random Batches. This figure illustrates the expression variability of the top 10 most differentially expressed genes (DEGs) in selected BMAC cell types across randomly assigned batches of 20 samples. Each dot represents a gene (y-axis) within a specific cell type (x-axis). Dot size reflects the percentage of cells expressing the gene, and color intensity denotes the average expression level, with red indicating higher expression and blue indicating lower expression. (A-D) The selected BMAC cell types display expression diversity for the most differentially expressed 10 genes, highlighting variability in gene expression patterns across sample batches. Cell-to-Cell Communication Analysis in BMAC and SVF Cell-to-cell communication analysis was conducted to ask whether there is potential for interactions between mesenchymal stromal cells (MSCs) and other cell types to differ between responders and non-responders. Using the CellChat tool, which infers likely receptor-ligand interaction networks from scRNA-seq data, a rich network of interactions was observed, with MSCs potentially capable of establishing communicative links with every cell type analyzed, and significant autocrine interactions noted in every cell type except macrophages ( Figure S3A ). Focusing on the top ten differentially expressed ligand-receptor interactions ( Figure S3B ), it was found that MSCs were pivotal predicted communicators in responder-derived BMAC cells. The relative information flow in Figure S3C suggests that MSCs in responder samples exhibited interactions enriched for ligands such as THY1, CD40, and IGF1, which are associated with immune modulation and tissue repair. By contrast, non-responder samples seem to be enriched for ligands such as IGFBP and CSF that are more typically linked to inflammatory and cellular stress responses. A similar analysis was performed on the SVF single cell RNAseq dataset, contrasting predicted cell communication differences between non-responders as the reference group and responders as the target condition. This analysis revealed ASC, MSC and M2 macrophages as the most prominent signal sender cell types in SVF ( Figure 7A ). The identities of the ASC and MSC populations, expression of gene markers shown in Figure 7B was used to separate ASCs and MSCs based on distinct expression profiles. MSCs exhibit higher expression of THY1 (CD90) and ENG , while ASCs are characterized by reduced CD34, HLA-DRA , and HLA-DRB1 expression. Focusing on specific ligand-receptor interactions and target gene predictions, several ligands enriched in the responder group were identified, both in terms of the number of interactions and their relative strengths. Dot plots in Figure 7C illustrate how the expression of various ligands differ between the three cell types. Notably, FGF7 is significantly upregulated in MSCs but is downregulated in ASCs, and conversely, IL15 is downregulated in MSCs and upregulated in ASCs. Other significant ligands enriched in MSCs include COL5A2, ANGPT1, POSTN, LAMA4 , and S100A4 , while ligands such as F8, NRG2 , and VWF are upregulated in ASCs. Additionally, the heatmap highlights the differential expression of genes predicted by NicheNet to be targets of the ligands, including CCL3, CCL4, IL2RB, STAT4, TRAF1, SESN1 , and NFKB1 . Download figure Open in new tab Figure 7: Interaction of Ligand Activity with Gene Expression and Predicted Regulatory Potential in SVF. A . Overall Interactions in SVF celltypes, this chord diagram illustrates the comprehensive landscape of cell-to-cell communication within SVF tissue. Each node represents a different cell type, with the connecting edges visualizing the complex web of potential signaling pathways. The diversity of edge thicknesses reflects the varying strengths of these interactions, providing insight into the intricate signaling network that orchestrates cellular functions and responses in SVF. B . Violin plot showing the ASC and MSC specific marker genes to help differentiate adipose stromal cells and mesenchymal stromal cells in the SVF tissue niche. C . The dot plot represents ligand activity across various prioritized ligands. The size of each dot corresponds to the percentage of cells expressing a given ligand within the sender population, with larger dots indicating a higher percentage. The color intensity reflects the average expression level of each ligand, with warmer colors (red-orange) representing higher expression levels in responders and cooler colors (yellow) representing lower levels in responders. Expanding on the analysis, a complementary investigation was conducted to examine the role of M2 macrophages as signal senders in SVF, adding another dimension to the predicted cell communication landscape ( Figure 7D ). The dot plot illustrates distinct ligand expression profiles in M2 macrophages compared to ASCs and MSCs, highlighting key differences in signaling activity. Notably, IFITM1 emerges as a significantly upregulated ligand in M2 macrophages, whereas its expression remains relatively lower in the other two cell types. Similarly, HBEGF, TGFB2, and FGF7 demonstrate increased expression in M2 macrophages, reinforcing their potential role in modulating the responder group’s microenvironment. In contrast, ligands such as CD86, BACE2, and VWF are more prominently expressed in ASCs and MSCs, suggesting differential contributions to the signaling network. This extended analysis underscores the nuanced cellular interactions within SVF and refines our understanding of how distinct cell populations contribute to the observed communication dynamics. Discussion Our study provides a comprehensive transcriptomic analysis of Bone Marrow Aspirate Concentrate (BMAC) and Stromal Vascular Fraction (SVF) samples in the context of osteoarthritis (OA) therapy. By analyzing cellular diversity and molecular expression patterns, our purpose was to identify key factors influencing treatment responses. We explored cell composition, gene expression variability, and cell communication networks, building on our previous work with a focus on treatment response variability. Our findings reveal marked study site-influenced variability in BMAC, with significant differences across five study sites in the United States affecting both cell composition and gene expression. This variability was more pronounced in BMAC than SVF, likely due to the heterogeneity incurred by the technique of needle aspiration from the bone marrow compared with more uniform sampling from subcutaneous fat. Operational variability may be unavoidable in real-world settings, but this finding nevertheless emphasizes the importance of testing and optimizing normalization techniques to enhance the likelihood that observed patterns reflect biological variability rather than procedural artifacts. We applied a range of normalization methods, including Harmony and rPCA, to adjust for site-related effects while preserving biological variability, showing how application of normalization techniques can potentially influence data interpretation. Overall, there was no significant enrichment of differential gene expression between responders and non-responders after adjusting for the multiple cell types and number of genes considered. Nevertheless, through scITD analysis, we were able to identify shared pathways of co-regulated genes across cell types in both BMAC and SVF. In BMAC, further implementation of a likelihood ratio test provided some evidence for key pathways that larger sample sizes may show are associated with treatment response. These included oxidative phosphorylation 30 , unfolded protein response 31 , and TNFα signaling 32 , suggestive of the initiation of a senescent 33 , pathophysiological environment 34 , related to cellular stress 35 and immune regulation 36 . Notably, the unfolded protein response, which is linked to endoplasmic reticulum stress, has been implicated in OA 37 , providing a potential mechanistic connection to disease progression. These findings suggest a heightened cellular stress response and metabolic dysregulation in non-responders, possibly contributing to the diminished therapeutic efficacy observed in these patients. In contrast, SVF samples showed distinct immune regulatory profiles 38 , 39 , with stromal cells exhibiting the potential to activate pathways related to innate immune responses and T cell activation 38 , while immune cells were enriched in cytokine signaling and antigen presentation pathways. These differences suggest that while BMAC cells may primarily modulate metabolic and inflammatory responses, SVF cells may have more potential to contribute to immune modulation. This dual role of metabolic dysregulation in BMAC and immune regulation 40 - 42 in SVF underscores the complexity of OA pathophysiology and highlights the need to tailor cell-based therapies based on these distinct molecular profiles. The pathway distinctions observed between responders and non-responders in BMAC and the immune-modulatory roles of SVF further emphasize the importance of understanding cellular contributions to therapeutic responses in OA 43 . The reason why the standard differential expression analysis did not reveal any statistically significant genes (after FDR correction) between responders and non-responders in either BMAC or SVF, whereas the LRT showed a strong pattern of enrichment, has to do with the combination of individual variability and study site confounding. The subset of individuals who failed to complete the study included many who dropped out due to non-response, explaining why their baseline profiles appear to be even more deviant from the healthy controls than the included non-responders. We expect that larger sample sizes will confirm the potential markers of non-response, and better address how to avoid study-site differences. Cell-cell communication analysis further underscored the potential therapeutic roles of MSCs within BMAC and SVF samples. In BMAC, MSCs were identified as major communicators, interacting with various immune cells. Notably, MSCs from responders showed enriched ligand interactions with THY1 44 , CD40 , and IGF1 45 , which are linked to immune modulation and tissue repair, whereas MSCs from non-responders exhibited elevated HLA signaling 46 , indicative of an inflammatory and possibly an autoimmune profile. This suggests that MSC-mediated signaling may be more coordinated in some responders, potentially contributing to treatment efficacy through enhanced immune cell interactions. In SVF, both MSCs and adipose-derived stromal cells (ASCs) demonstrated active communication networks, particularly in responders. Ligands such as IL15, FGF7 , and TNFSF12 showed varying expression patterns 47 - 49 , with FGF7 upregulated in MSCs and downregulated in ASCs 50 , indicating nuanced differences in cellular interactions that could affect immune responses and tissue repair in OA 51 , 52 . All of these cell communication inferences are subject to inter-individual variability. These findings underscore a major limitation of transcriptomic profiling in the context of therapeutic responses, which is that cell therapy outcomes in OA are subject to both patient-specific factors and site-dependent variability. While genomic profiling holds potential for understanding treatment mechanisms, our results indicate that it may be premature to guide OA therapy based solely on molecular profiling, given the high degree of variability. The suggestive pathways identified in this study provide a foundation for future research to clarify these molecular mechanisms and explore their role in patient response, paving the way for more targeted cell-based therapies in OA. Overall, this study highlights the importance of accounting for site variability in BMAC samples, the distinct roles of MSCs in BMAC and SVF, and the potential of integrating scRNA-seq and cell communication analyses to refine therapeutic strategies. Our findings contribute to a growing understanding of OA’s cellular dynamics and suggest that personalized cell therapy may require a more nuanced approach, considering both molecular signatures and sample-specific factors to improve treatment outcomes. Data Availability All data produced in the present study will be available online once published in a peer-reviewed journal. Supplemental Information Supplemental Information can be found online at [Link will be provided after publishing with peer-reviewed journal] Data and Code availability The datasets used in this manuscript are available with a reviewer link to GSE292453. Will be made public upon publishing this paper. Author Contributions P.C. performed scRNAseq experiments and analyzed and interpreted the data with guidance from G.G. H.S. and L.K. assisted with sample processing and quality control. S.B., K.M. and H.D. managed the MILES trial and conceived the overall experimental design. C.Y., H.D., and K.R. also contributed to data interpretation. All authors reviewed and approved the final manuscript which was initially drafted by P.C. and G.G. Conflicts of interest The authors declare no competing interests. Acknowledgements We wish to acknowledge the Cellular Analysis and Cytometry Core and the Molecular Evolution Core facilities at the Parker H. Petit Institute for Bioengineering and Bioscience at the 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 Marcus Foundation to Emory University and the 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 Georgia Tech Foundation and the Georgia Research Alliance. References 1. ↵ Yao , Q. , Wu , X. , Tao , C. , Gong , W. , Chen , M. , Qu , M. , Zhong , Y. , He , T. , Chen , S. , and Xiao , G. ( 2023 ). Osteoarthritis: pathogenic signaling pathways and therapeutic targets . Signal Transduct Target Ther 8 , 56 . doi: 10.1038/s41392-023-01330-w . OpenUrl CrossRef 2. ↵ Chen , T. , Weng , W. , Liu , Y. , Aspera-Werz , R.H. , Nussler , A.K. , and Xu , J. ( 2021 ). Update on Novel Non-Operative Treatment for Osteoarthritis: Current Status and Future Trends . Front Pharmacol 12 , 755230 . doi: 10.3389/fphar.2021.755230 . OpenUrl CrossRef PubMed 3. ↵ Kim , G.B. , Seo , M.S. , Park , W.T. , and Lee , G.W. ( 2020 ). Bone Marrow Aspirate Concentrate: Its Uses in Osteoarthritis . Int J Mol Sci 21 . doi: 10.3390/ijms21093224 . OpenUrl CrossRef 4. ↵ Goncharov , E.N. , Koval , O.A. , Nikolaevich Bezuglov , E. , Encarnacion Ramirez , M.J. , Engelgard , M. , Igorevich , E.I. , Saporiti , A. , Valentinovich Kotenko , K. , and Montemurro , N. ( 2023 ). Stromal Vascular Fraction Therapy for Knee Osteoarthritis: A Systematic Review . Medicina (Kaunas) 59 . doi: 10.3390/medicina59122090 . OpenUrl CrossRef 5. ↵ Michalek , J. , Vrablikova , A. , Darinskas , A. , Lukac , L. , Prucha , J. , Skopalik , J. , Travnik , J. , Cibulka , M. , and Dudasova , Z. ( 2019 ). Stromal vascular fraction cell therapy for osteoarthritis in elderly: Multicenter case-control study . J Clin Orthop Trauma 10 , 76 – 80 . doi: 10.1016/j.jcot.2018.11.010 . OpenUrl CrossRef PubMed 6. ↵ Silverwood , V. , Blagojevic-Bucknall , M. , Jinks , C. , Jordan , J.L. , Protheroe , J. , and Jordan , K.P. ( 2015 ). Current evidence on risk factors for knee osteoarthritis in older adults: a systematic review and meta-analysis . Osteoarthritis Cartilage 23 , 507 – 515 . doi: 10.1016/j.joca.2014.11.019 . OpenUrl CrossRef PubMed 7. ↵ Mautner , K. , Gottschalk , M. , Boden , S.D. , Akard , A. , Bae , W.C. , Black , L. , Boggess , B. , Chatterjee , P. , Chung , C.B. , Easley , K.A. , Gibson , G. , et al. ( 2024 ). Author Correction: Cell-based versus corticosteroid injections for knee pain in osteoarthritis: a randomized phase 3 trial . Nat Med 30 , 905 . doi: 10.1038/s41591-023-02776-9 . OpenUrl CrossRef 8. ↵ Chatterjee , P. , Stevens , H.Y. , Kippner , L.E. , Bowles-Welch , A.C. , Drissi , H. , Mautner , K. , Yeago , C. , Gibson , G. , and Roy , K. ( 2024 ). Single-cell transcriptome and crosstalk analysis reveals immune alterations and key pathways in the bone marrow of knee OA patients . iScience 27 , 110827 . doi: 10.1016/j.isci.2024.110827 . OpenUrl CrossRef PubMed 9. ↵ Sebastian , A. , Hum , N.R. , McCool , J.L. , Wilson , S.P. , Murugesh , D.K. , Martin , K.A. , Rios-Arce , N.D. , Amiri , B. , Christiansen , B.A. , and Loots , G.G. ( 2022 ). Single-cell RNA-Seq reveals changes in immune landscape in post-traumatic osteoarthritis . Front Immunol 13 , 938075 . doi: 10.3389/fimmu.2022.938075 . OpenUrl CrossRef PubMed 10. ↵ Wang , Y. , Wluka , A.E. , Jones , G. , Ding , C. , and Cicuttini , F.M. ( 2012 ). Use magnetic resonance imaging to assess articular cartilage . Ther Adv Musculoskelet Dis 4 , 77 – 97 . doi: 10.1177/1759720X11431005 . OpenUrl CrossRef PubMed 11. ↵ Kang , M. , Ragan , B.G. , and Park , J.H. ( 2008 ). Issues in outcomes research: an overview of randomization techniques for clinical trials . J Athl Train 43 , 215 – 221 . doi: 10.4085/1062-6050-43.2.215 . OpenUrl CrossRef PubMed Web of Science 12. ↵ Tran , H.T.N. , Ang , K.S. , Chevrier , M. , Zhang , X. , Lee , N.Y.S. , Goh , M. , and Chen , J. ( 2020 ). A benchmark of batch-effect correction methods for single-cell RNA sequencing data . Genome Biol 21 , 12 . doi: 10.1186/s13059-019-1850-9 . OpenUrl CrossRef PubMed 13. ↵ Hao , Y. , Hao , S. , Andersen-Nissen , E. , Mauck , W.M. , 3rd . , Zheng , S. , Butler , A. , Lee , M.J. , Wilk , A.J. , Darby , C. , Zager , M. , Hoffman , P. , et al. ( 2021 ). Integrated analysis of multimodal single-cell data . Cell 184 , 3573 – 3587 e3529 . doi: 10.1016/j.cell.2021.04.048 . OpenUrl CrossRef PubMed 14. ↵ Korsunsky , I. , Millard , N. , Fan , J. , Slowikowski , K. , Zhang , F. , Wei , K. , Baglaenko , Y. , Brenner , M. , Loh , P.R. , and Raychaudhuri , S. ( 2019 ). Fast, sensitive and accurate integration of single-cell data with Harmony . Nat Methods 16 , 1289 – 1296 . doi: 10.1038/s41592-019-0619-0 . OpenUrl CrossRef PubMed 15. ↵ Li , M. , Zhang , X. , Ang , K.S. , Ling , J. , Sethi , R. , Lee , N.Y.S. , Ginhoux , F. , and Chen , J. ( 2022 ). DISCO: a database of Deeply Integrated human Single-Cell Omics data . Nucleic Acids Res 50 , D596 – D602 . doi: 10.1093/nar/gkab1020 . OpenUrl CrossRef 16. ↵ Becht , E. , McInnes , L. , Healy , J. , Dutertre , C.A. , Kwok , I.W.H. , Ng , L.G. , Ginhoux , F. , and Newell , E.W. ( 2018 ). Dimensionality reduction for visualizing single-cell data using UMAP . Nat Biotechnol . doi: 10.1038/nbt.4314 . OpenUrl CrossRef PubMed 17. ↵ Phipson , B. , Sim , C.B. , Porrello , E.R. , Hewitt , A.W. , Powell , J. , and Oshlack , A. ( 2022 ). propeller: testing for differences in cell type proportions in single cell data . Bioinformatics 38 , 4720 – 4726 . doi: 10.1093/bioinformatics/btac582 . OpenUrl CrossRef PubMed 18. ↵ Murphy , A.E. , and Skene , N.G. ( 2022 ). A balanced measure shows superior performance of pseudobulk methods in single-cell RNA-sequencing analysis . Nat Commun 13 , 7851 . doi: 10.1038/s41467-022-35519-4 . OpenUrl CrossRef PubMed 19. ↵ Zhang , Y. , Parmigiani , G. , and Johnson , W.E. ( 2020 ). ComBat-seq: batch effect adjustment for RNA-seq count data . NAR Genom Bioinform 2 , lqaa078 . doi: 10.1093/nargab/lqaa078 . OpenUrl CrossRef PubMed 20. ↵ Love , M.I. , Huber , W. , and Anders , S. ( 2014 ). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 . Genome Biol 15 , 550 . doi: 10.1186/s13059-014-0550-8 . OpenUrl CrossRef PubMed 21. ↵ Perneger , T.V. ( 2021 ). How to use likelihood ratios to interpret evidence from randomized trials . J Clin Epidemiol 136 , 235 – 242 . doi: 10.1016/j.jclinepi.2021.04.010 . OpenUrl CrossRef PubMed 22. ↵ Hoffman , G.E. , Lee , D. , Bendl , J. , Fnu , P. , Hong , A. , Casey , C. , Alvia , M. , Shao , Z. , Argyriou , S. , Therrien , K. , Venkatesh , S. , et al. ( 2023 ). Efficient differential expression analysis of large-scale single cell transcriptomics data using dreamlet . Res Sq . doi: 10.21203/rs.3.rs-2705625/v1 . OpenUrl CrossRef 23. ↵ Zimmerman , K.D. , Evans , C. , and Langefeld , C.D. ( 2022 ). Reply to: A balanced measure shows superior performance of pseudobulk methods in single-cell RNA-sequencing analysis . Nat Commun 13 , 7852 . doi: 10.1038/s41467-022-35520-x . OpenUrl CrossRef PubMed 24. ↵ Chen , X. , Doerge , R.W. , and Sarkar , S.K. ( 2020 ). A weighted FDR procedure under discrete and heterogeneous null distributions . Biom J 62 , 1544 – 1563 . doi: 10.1002/bimj.201900216 . OpenUrl CrossRef PubMed 25. ↵ Jin , S. , Guerrero-Juarez , C.F. , Zhang , L. , Chang , I. , Ramos , R. , Kuan , C.H. , Myung , P. , Plikus , M.V. , and Nie , Q. ( 2021 ). Inference and analysis of cell-cell communication using CellChat . Nat Commun 12 , 1088 . doi: 10.1038/s41467-021-21246-9 . OpenUrl CrossRef PubMed 26. ↵ Browaeys , R. , Saelens , W. , and Saeys , Y. ( 2020 ). NicheNet: modeling intercellular communication by linking ligands to target genes . Nat Methods 17 , 159 – 162 . doi: 10.1038/s41592-019-0667-5 . OpenUrl CrossRef PubMed 27. ↵ Subramanian , A. , Tamayo , P. , Mootha , V.K. , Mukherjee , S. , Ebert , B.L. , Gillette , M.A. , Paulovich , A. , Pomeroy , S.L. , Golub , T.R. , Lander , E.S. , and Mesirov , J.P. ( 2005 ). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles . Proc Natl Acad Sci U S A 102 , 15545 – 15550 . doi: 10.1073/pnas.0506580102 . OpenUrl Abstract / FREE Full Text 28. ↵ Ma , Y. , Sun , S. , Shang , X. , Keller , E.T. , Chen , M. , and Zhou , X. ( 2020 ). Integrative differential expression and gene set enrichment analysis using summary statistics for scRNA-seq studies . Nat Commun 11 , 1585 . doi: 10.1038/s41467-020-15298-6 . OpenUrl CrossRef PubMed 29. ↵ Andreatta , M. , and Carmona , S.J. ( 2021 ). UCell: Robust and scalable single-cell gene signature scoring . Comput Struct Biotechnol J 19 , 3796 – 3798 . doi: 10.1016/j.csbj.2021.06.043 . OpenUrl CrossRef PubMed 30. ↵ Han , J. , Kim , Y.H. , and Han , S. ( 2025 ). Increased oxidative phosphorylation through pyruvate dehydrogenase kinase 2 deficiency ameliorates cartilage degradation in mice with surgically induced osteoarthritis . Exp Mol Med 57 , 390 – 401 . doi: 10.1038/s12276-025-01400-9 . OpenUrl CrossRef PubMed 31. ↵ Li , Y.H. , Tardif , G. , Hum , D. , Kapoor , M. , Fahmi , H. , Pelletier , J.P. , and Martel-Pelletier , J. ( 2016 ). The unfolded protein response genes in human osteoarthritic chondrocytes: PERK emerges as a potential therapeutic target . Arthritis Res Ther 18 , 172 . doi: 10.1186/s13075-016-1070-6 . OpenUrl CrossRef PubMed 32. ↵ Cavallo , C. , Boffa , A. , Andriolo , L. , Silva , S. , Grigolo , B. , Zaffagnini , S. , and Filardo , G. ( 2021 ). Bone marrow concentrate injections for the treatment of osteoarthritis: evidence from preclinical findings to the clinical application . Int Orthop 45 , 525 – 538 . doi: 10.1007/s00264-020-04703-w . OpenUrl CrossRef PubMed 33. ↵ Jeon , O.H. , David , N. , Campisi , J. , and Elisseeff , J.H. ( 2018 ). Senescent cells and osteoarthritis: a painful connection . J Clin Invest 128 , 1229 – 1237 . doi: 10.1172/JCI95147 . OpenUrl CrossRef PubMed 34. ↵ Yunus , M.H.M. , Nordin , A. , and Kamal , H. ( 2020 ). Pathophysiological Perspective of Osteoarthritis . Medicina (Kaunas) 56 . doi: 10.3390/medicina56110614 . OpenUrl CrossRef 35. ↵ Chen , Y. , Cheng , R.J. , Wu , Y. , Huang , D. , Li , Y. , and Liu , Y. ( 2023 ). Advances in Stem Cell-Based Therapies in the Treatment of Osteoarthritis . Int J Mol Sci 25 . doi: 10.3390/ijms25010394 . OpenUrl CrossRef 36. ↵ Pabinger , C. , Lothaller , H. , and Kobinia , G.S. ( 2024 ). Intra-articular injection of bone marrow aspirate concentrate (mesenchymal stem cells) in KL grade III and IV knee osteoarthritis: 4 year results of 37 knees . Sci Rep 14 , 2665 . doi: 10.1038/s41598-024-51410-2 . OpenUrl CrossRef PubMed 37. ↵ Soh , L.J. , Lee , S.Y. , Roebuck , M.M. , and Wong , P.F. ( 2024 ). Unravelling the interplay between ER stress, UPR and the cGAS-STING pathway: Implications for osteoarthritis pathogenesis and treatment strategy . Life Sci 357 , 123112 . doi: 10.1016/j.lfs.2024.123112 . OpenUrl CrossRef PubMed 38. ↵ Goncharov , E.N. , Koval , O.A. , Nikolaevich Bezuglov , E. , Engelgard , M. , Igorevich , E.I. , Velentinovich Kotenko , K. , Encarnacion Ramirez , M.J. , and Montemurro , N. ( 2024 ). Comparative Analysis of Stromal Vascular Fraction and Alternative Mechanisms in Bone Fracture Stimulation to Bridge the Gap between Nature and Technological Advancement: A Systematic Review . Biomedicines 12 . doi: 10.3390/biomedicines12020342 . OpenUrl CrossRef 39. ↵ Onoi , Y. , Matsumoto , T. , Anjiki , K. , Hayashi , S. , Nakano , N. , Kuroda , Y. , Tsubosaka , M. , Kamenaga , T. , Ikuta , K. , Tachibana , S. , Suda , Y. , et al. ( 2024 ). Human uncultured adipose-derived stromal vascular fraction shows therapeutic potential against osteoarthritis in immunodeficient rats via direct effects of transplanted M2 macrophages . Stem Cell Res Ther 15 , 325 . doi: 10.1186/s13287-024-03946-3 . OpenUrl CrossRef PubMed 40. ↵ Liu , T. , Zhang , L. , Joo , D. , and Sun , S.C. ( 2017 ). NF-kappaB signaling in inflammation . Signal Transduct Target Ther 2 , 17023 -. doi: 10.1038/sigtrans.2017.23 . OpenUrl CrossRef PubMed 41. Furman , B.D. , Kent , C.L. , Huebner , J.L. , Kraus , V.B. , McNulty , A.L. , Guilak , F. , and Olson , S.A. ( 2018 ). CXCL10 is upregulated in synovium and cartilage following articular fracture . J Orthop Res 36 , 1220 – 1227 . doi: 10.1002/jor.23735 . OpenUrl CrossRef PubMed 42. ↵ Hari , P. , Millar , F.R. , Tarrats , N. , Birch , J. , Quintanilla , A. , Rink , C.J. , Fernandez-Duran , I. , Muir , M. , Finch , A.J. , Brunton , V.G. , Passos , J.F. , et al. ( 2019 ). The innate immune sensor Toll-like receptor 2 controls the senescence-associated secretory phenotype . Sci Adv 5 , eaaw0254 . doi: 10.1126/sciadv.aaw0254 . OpenUrl FREE Full Text 43. ↵ Lee , Y.T. , Mohd Yunus , M.H. , Yazid , M.D. , and Ugusman , A. ( 2024 ). Unraveling the path to osteoarthritis management: targeting chondrocyte apoptosis for therapeutic intervention . Front Cell Dev Biol 12 , 1347126 . doi: 10.3389/fcell.2024.1347126 . OpenUrl CrossRef PubMed 44. ↵ Paine , A. , Woeller , C.F. , Zhang , H. , de la Luz Garcia-Hernandez , M. , Huertas , N. , Xing , L. , Phipps , R.P. , and Ritchlin , C.T. ( 2018 ). Thy1 is a positive regulator of osteoblast differentiation and modulates bone homeostasis in obese mice . FASEB J 32 , 3174 – 3183 . doi: 10.1096/fj.201701379R . OpenUrl CrossRef 45. ↵ Wu , H. , Peng , Z. , Xu , Y. , Sheng , Z. , Liu , Y. , Liao , Y. , Wang , Y. , Wen , Y. , Yi , J. , Xie , C. , Chen , X. , et al. ( 2022 ). Engineered adipose-derived stem cells with IGF-1-modified mRNA ameliorates osteoarthritis development . Stem Cell Res Ther 13 , 19 . doi: 10.1186/s13287-021-02695-x . OpenUrl CrossRef PubMed 46. ↵ Wang , Y. , Huang , J. , Gong , L. , Yu , D. , An , C. , Bunpetch , V. , Dai , J. , Huang , H. , Zou , X. , Ouyang , H. , and Liu , H. ( 2019 ). The Plasticity of Mesenchymal Stem Cells in Regulating Surface HLA-I . iScience 15 , 66 – 78 . doi: 10.1016/j.isci.2019.04.011 . OpenUrl CrossRef PubMed 47. ↵ Gonzalez-Cubero , E. , Gonzalez-Fernandez , M.L. , Esteban-Blanco , M. , Perez-Castrillo , S. , Perez-Fernandez , E. , Navasa , N. , Aransay , A.M. , Anguita , J. , and Villar-Suarez , V. ( 2024 ). The Therapeutic Potential of Adipose-Derived Mesenchymal Stem Cell Secretome in Osteoarthritis: A Comprehensive Study . Int J Mol Sci 25 . doi: 10.3390/ijms252011287 . OpenUrl CrossRef 48. Zhang , X. , Pu , X. , Pi , C. , and Xie , J. ( 2023 ). The role of fibroblast growth factor 7 in cartilage development and diseases . Life Sci 326 , 121804 . doi: 10.1016/j.lfs.2023.121804 . OpenUrl CrossRef PubMed 49. ↵ Matellan , C. , Kennedy , C. , Santiago-Vela , M.I. , Hochegger , J. , Ni Chathail , M.B. , Wu , A. , Shannon , C. , Roche , H.M. , Aceves , S.S. , Godson , C. , and Manresa , M.C. ( 2024 ). The TNFSF12/TWEAK Modulates Colonic Inflammatory Fibroblast Differentiation and Promotes Fibroblast-Monocyte Interactions . J Immunol 212 , 1958 – 1970 . doi: 10.4049/jimmunol.2300762 . OpenUrl CrossRef PubMed 50. ↵ Shahin , H. , Belcastro , L. , Das , J. , Perdiki Grigoriadi , M. , Saager , R.B. , Steinvall , I. , Sjoberg , F. , Olofsson , P. , Elmasry , M. , and El-Serafi , A.T. ( 2024 ). MicroRNA-155 mediates multiple gene regulations pertinent to the role of human adipose-derived mesenchymal stem cells in skin regeneration . Front Bioeng Biotechnol 12 , 1328504 . doi: 10.3389/fbioe.2024.1328504 . OpenUrl CrossRef PubMed 51. ↵ Wei , J. , and Besner , G.E. ( 2015 ). M1 to M2 macrophage polarization in heparin-binding epidermal growth factor-like growth factor therapy for necrotizing enterocolitis . J Surg Res 197 , 126 – 138 . doi: 10.1016/j.jss.2015.03.023 . OpenUrl CrossRef PubMed 52. ↵ Becker , M. , De Bastiani , M.A. , Parisi , M.M. , Guma , F.T. , Markoski , M.M. , Castro , M.A. , Kaplan , M.H. , Barbe-Tuana , F.M. , and Klamt , F. ( 2015 ). Integrated Transcriptomics Establish Macrophage Polarization Signatures and have Potential Applications for Clinical Health and Disease . Sci Rep 5 , 13351 . doi: 10.1038/srep13351 . OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted April 16, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about medRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Single-Cell Transcriptomics of Multi-Site Cell Therapy in Osteoarthritis: Tissue-Specific Traits and Treatment Correlations Message Subject (Your Name) has forwarded a page to you from medRxiv Message Body (Your Name) thought you would like to see this page from the medRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Single-Cell Transcriptomics of Multi-Site Cell Therapy in Osteoarthritis: Tissue-Specific Traits and Treatment Correlations Paramita Chatterjee , Hazel Y. Stevens , Linda E. Kippner , Carolyn Yeago , Hicham Drissi , Kenneth Mautner , Scott Boden , Greg Gibson , Krishnendu Roy medRxiv 2025.04.14.25325743; doi: https://doi.org/10.1101/2025.04.14.25325743 Share This Article: Copy Citation Tools Single-Cell Transcriptomics of Multi-Site Cell Therapy in Osteoarthritis: Tissue-Specific Traits and Treatment Correlations Paramita Chatterjee , Hazel Y. Stevens , Linda E. Kippner , Carolyn Yeago , Hicham Drissi , Kenneth Mautner , Scott Boden , Greg Gibson , Krishnendu Roy medRxiv 2025.04.14.25325743; doi: https://doi.org/10.1101/2025.04.14.25325743 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Genetic and Genomic Medicine Subject Areas All Articles Addiction Medicine (568) Allergy and Immunology (863) Anesthesia (300) Cardiovascular Medicine (4435) Dentistry and Oral Medicine (444) Dermatology (382) Emergency Medicine (608) Endocrinology (including Diabetes Mellitus and Metabolic Disease) (1509) Epidemiology (15229) Forensic Medicine (30) Gastroenterology (1124) Genetic and Genomic Medicine (6600) Geriatric Medicine (668) Health Economics (997) Health Informatics (4536) Health Policy (1368) Health Systems and Quality Improvement (1613) Hematology (541) HIV/AIDS (1264) Infectious Diseases (except HIV/AIDS) (15916) Intensive Care and Critical Care Medicine (1103) Medical Education (623) Medical Ethics (146) Nephrology (667) Neurology (6599) Nursing (346) Nutrition (998) Obstetrics and Gynecology (1144) Occupational and Environmental Health (957) Oncology (3332) Ophthalmology (974) Orthopedics (369) Otolaryngology (420) Pain Medicine (436) Palliative Medicine (130) Pathology (663) Pediatrics (1693) Pharmacology and Therapeutics (691) Primary Care Research (711) Psychiatry and Clinical Psychology (5447) Public and Global Health (9232) Radiology and Imaging (2198) Rehabilitation Medicine and Physical Therapy (1370) Respiratory Medicine (1196) Rheumatology (593) Sexual and Reproductive Health (712) Sports Medicine (530) Surgery (712) Toxicology (99) Transplantation (289) Urology (265) (function(){function c(){var b=a.contentDocument||a.contentWindow.document;if(b){var d=b.createElement('script');d.innerHTML="window.__CF$cv$params={r:'a00ad7b1ab790db4',t:'MTc3OTYxMDQyMg=='};var a=document.createElement('script');a.src='/cdn-cgi/challenge-platform/scripts/jsd/main.js';document.getElementsByTagName('head')[0].appendChild(a);";b.getElementsByTagName('head')[0].appendChild(d)}}if(document.body){var a=document.createElement('iframe');a.height=1;a.width=1;a.style.position='absolute';a.style.top=0;a.style.left=0;a.style.border='none';a.style.visibility='hidden';document.body.appendChild(a);if('loading'!==document.readyState)c();else if(window.addEventListener)document.addEventListener('DOMContentLoaded',c);else{var e=document.onreadystatechange||function(){};document.onreadystatechange=function(b){e(b);'loading'!==document.readyState&&(document.onreadystatechange=e,c())}}}})();