Transcriptional and Histochemical Signatures of Bone Marrow Mononuclear Cell-Mediated Resolution of Synovitis

Osteoarthritis (OA) may result from impaired ability of synovial macrophages to resolve joint inflammation. Increasing macrophage counts in inflamed joints through injection with bone marrow mononuclear cells (BMNC) induces lasting resolution of synovial inflammation. To uncover mechanisms by which BMNC may affect resolution, in this study, differential transcriptional signatures of BMNC in response to normal (SF) and inflamed synovial fluid (ISF) were analyzed. We demonstrate the temporal behavior of co-expressed gene networks associated with traits from related in vivo and in vitro studies. We also identified activated and inhibited signaling pathways and upstream regulators, further determining their protein expression in the synovium of inflamed joints treated with BMNC or DPBS controls. BMNC responded to ISF with an early pro-inflammatory response characterized by a short spike in the expression of a NF-ƙB- and mitogen-related gene network. This response was associated with sustained increased expression of two gene networks comprising known drivers of resolution (IL-10, IGF-1, PPARG, isoprenoid biosynthesis). These networks were common to SF and ISF, but more highly expressed in ISF. Most highly activated pathways in ISF included the mevalonate pathway and PPAR-γ signaling, with pro-resolving functional annotations that improve mitochondrial metabolism and deactivate NF-ƙB signaling. Lower expression of mevalonate kinase and phospho-PPARγ in synovium from inflamed joints treated with BMNC, and equivalent IL-1β staining between BMNC- and DPBS-treated joints, associates with accomplished resolution in BMNC-treated joints and emphasize the intricate balance of pro- and anti-inflammatory mechanisms required for resolution. Combined, our data suggest that BMNC-mediated resolution is characterized by constitutively expressed homeostatic mechanisms, whose expression are enhanced following inflammatory stimulus. These mechanisms translate into macrophage proliferation optimizing their capacity to counteract inflammatory damage and improving their general and mitochondrial metabolism to endure oxidative stress while driving tissue repair. Such effect is largely achieved through the synthesis of several lipids that mediate recovery of homeostasis. Our study reveals candidate mechanisms by which BMNC provide lasting improvement in patients with OA and suggests further investigation on the effects of PPAR-γ signaling enhancement for the treatment of arthritic conditions.

Osteoarthritis (OA) may result from impaired ability of synovial macrophages to resolve joint inflammation. Increasing macrophage counts in inflamed joints through injection with bone marrow mononuclear cells (BMNC) induces lasting resolution of synovial inflammation. To uncover mechanisms by which BMNC may affect resolution, in this study, differential transcriptional signatures of BMNC in response to normal (SF) and inflamed synovial fluid (ISF) were analyzed. We demonstrate the temporal behavior of coexpressed gene networks associated with traits from related in vivo and in vitro studies. We also identified activated and inhibited signaling pathways and upstream regulators, further determining their protein expression in the synovium of inflamed joints treated with BMNC or DPBS controls. BMNC responded to ISF with an early pro-inflammatory response characterized by a short spike in the expression of a NF-ƙB-and mitogenrelated gene network. This response was associated with sustained increased expression of two gene networks comprising known drivers of resolution (IL-10, IGF-1, PPARG, isoprenoid biosynthesis). These networks were common to SF and ISF, but more highly expressed in ISF. Most highly activated pathways in ISF included the mevalonate pathway and PPAR-g signaling, with pro-resolving functional annotations that improve mitochondrial metabolism and deactivate NF-ƙB signaling. Lower expression of mevalonate kinase and phospho-PPARg in synovium from inflamed joints treated with BMNC, and equivalent IL-1b staining between BMNC-and DPBS-treated joints, associates with accomplished resolution in BMNC-treated joints and emphasize the intricate balance of pro-and anti-inflammatory mechanisms required for resolution. Combined, our data suggest that BMNC-mediated resolution is characterized by constitutively expressed homeostatic mechanisms, whose expression are enhanced following inflammatory stimulus. These mechanisms translate into macrophage proliferation optimizing their capacity to counteract inflammatory damage and improving their general and mitochondrial metabolism to endure oxidative stress while driving tissue INTRODUCTION Osteoarthritis (OA) is a common and debilitating condition that similarly affects horses and people (1,2). Because chronic synovial inflammation is a hallmark of OA and often the single driver of related degenerative changes (3)(4)(5)(6)(7), the use of antiinflammatory drugs (steroidal and non-steroidal) has been a logical and long-accepted approach for the treatment of many arthritic conditions (8,9). However, acute inflammation is not simply a clinical sign to alleviate. Acute inflammation is a critical event in promoting tissue repair and setting the stage for endogenous resolution of the inflammatory process and recovery of homeostasis (10). Importantly, anti-inflammatory and pro-resolving effects are not the same, and resolution is not merely the passive termination of the inflammatory process. Antiinflammation is based on inhibiting key pro-inflammatory mediators, such as chemokine and cytokine production and leukocyte extravasation to the site of injury. Resolution is an active process driven primarily by macrophages and their derived cytokines and lipid mediators, which shift the phlogistic phase of inflammation into a non-phlogistic process that culminates with tissue repair and recovery of homeostasis (11,12). Most importantly, the recruitment of macrophages and the production of pro-resolving mediators is triggered by enzymes synthesized during the acute inflammatory process (13). Macrophages play such a fundamental role in resolving inflammation and promoting tissue repair that impaired macrophage chemotaxis and/or macrophage depletion results in inefficient healing or chronic inflammation (14)(15)(16). Blocking acute inflammation with anti-inflammatory medications interferes, at least to some degree, with macrophage recruitment and the pro-resolving response, and often prevents effective resolution and recovery of homeostasis (12,13). Targeted therapies for chronic joint inflammation should therefore have pro-resolving properties, which precisely combine pro-and antiinflammatory mechanisms (12).
Synovial macrophages are the central drivers of the inflammatory response in osteoarthritic synovium (17,18). In fact, synovial macrophage activation is directly related to disease activity, severity, and pain in OA-affected patients (19). However, this relationship is not causative. Synovial macrophages are also essential keepers of synovial homeostasis through phagocytic clearance and secretion of anti-inflammatory and pro-resolving cytokines, chemokines, enzymes, and growth factors (20)(21)(22)(23). Following injury, synovial macrophages proliferate to form a protective immunological barrier in the synovial lining for intra-articular structures (24,25). When regulatory functions are overwhelmed by the amount of damage, synovial macrophages upregulate inflammation, signaling to monocytes and other leukocytes (e.g., neutrophils and lymphocytes) to help counteract the increased demands for tissue repair and restore homeostasis (17,26). During the progression of OA, the recruitment of myeloid monocytes into joints seems to be impaired (27), which combined with continuous joint damage, overwhelms the proresolving mechanisms of synovial macrophages, leading to degeneration (17,27,28).
The mononuclear cell fraction of bone marrow aspirates (bone marrow mononuclear cells -BMNC) is a rich source of proresolving macrophages that have been used therapeutically to improve tissue repair and inflammation resolution (29)(30)(31)(32)(33)(34)(35)(36)(37)(38)(39). Macrophages within BMNC are the main drivers of such effect, which's documented pro-resolving functions include increased production of IL-10 (30,33), diverse types of prostaglandins and specialized lipid mediators (13,40,41). The production of these molecules induce decreased production of IL-6 (33), increased phagocytic clearance of debris and apoptotic cells (efferocytosis) (30,39) and enhanced PPAR-gamma signaling (42)(43)(44)(45)(46)(47). Increasing the numbers of myeloid macrophages present in osteoarthritic knees, by injection of BMNC, restored joint homeostasis with long-lasting effects (48). Similarly, BMNC therapy increased counts of pro-resolving macrophages and induced marked resolution of joint inflammation in in vivo and in vitro models of equine synovitis (49,50). In these models, there was a coordinated spectrum of pro-inflammatory, pro-resolving and anti-inflammatory events, including increased IL-10, IGF-1, and PGE 2 production, and self-limiting IL-1 signaling (a and b). These events are all innately required for efficient synovial homeostasis and tissue repair and are commonly antagonized by therapeutic corticosteroids (30,(51)(52)(53). While these findings partially explain the durable effects of BMNC in the treatment of OA, little is known about BMNC-related mechanisms of resolution. Therefore, our purpose was to identify cellular mechanisms from BMNC driving joint homeostasis that could be used for developing targeted pro-resolving joint therapies and uncovering biomarkers of arthritis resolution. The aim of this study was to identify transcriptional signatures of BMNC leading to inflammation resolution using RNA-sequencing and relate these to the expression of key gene products in the synovial membrane. We hypothesized that gene networks linked to macrophage proliferation, negative regulation of inflammatory response, and to a lesser extent, NF-ƙB signaling, would be temporally upregulated in response to inflammation.

Study Design
Samples used in the current report were obtained from two previous studies using the same horses. These in vitro (50) and in vivo (49) studies were counterparts of a larger project assessing the effects of BMNC on joint inflammation resolution. Briefly, eight skeletally mature Thoroughbred horses (3-9 years old, median 5 years; 2 females and 6 castrated males) free of OA or systemic inflammation were used under IACUC approval and oversight. General and musculoskeletal health were confirmed by clinical, hematological and orthopedic evaluations. Following sternal bone marrow aspiration for BMNC isolation, synovitis was induced in both radiocarpal joints, as a way of producing more homogeneous inflammation and inflamed synovial fluid (ISF) than could be acquired from naturally occurring OA. Normal synovial fluid (SF) was collected from healthy middle carpal joints. BMNC from each horse were cultured independently (not pooled) in neat (100%) autologous SF or ISF and harvested at 0 (uncultured), 48 and 96 hours, and 6 and 10 days for RNA isolation. RNA-seq was used to identify transcriptional signatures of BMNC in response to acute joint inflammation. The transcriptome of BMNC was assessed over time within the same group (SF or ISF), as well as comparatively between groups at each time point (Figure 1). The expression of 9 potential upstream regulator genes identified following bioinformatical analysis was assessed by immunohistochemistry in the synovium of inflamed joints from the same horses, 6 days after treatment with BMNC or Dulbecco's phosphate buffered saline (DPBS). Four protein coded by genes identified in key activated pathways were also assessed by IHC.

BMNC Isolation, Induction of Synovitis, and Synovial Fluid Harvest
Bone marrow harvest and processing of BMNC, and induction of the synovitis model were performed as previously described in our related study (49). Briefly, BMNC were isolated by density gradient centrifugation. Synovitis was induced by intra-articular injection of 0.5 ng lipopolysaccharide (LPS) into each radiocarpal joint (49,54). At peak inflammation (8 hours following induction of the model), SF and ISF were collected using aseptic technique. Synovial fluid cytology (SF and ISF) was performed to confirm the health of normal joints and ensure LPS effectively induced synovitis. Synovial fluid was then centrifuged (5,000g; 20 min; 4°C) for cell depletion and the cell-free supernatant used as autologous growth medium. Parameters used to differentiate SF from ISF included quantification of cytokine as reported in our related study (50), total protein (<2.5 g/dL in SF and >4g/dL in ISF) and synovial fluid cytology (total nucleated cells/µL < 1,500 in SF and 130,000 in ISF; neutrophil count <10% in SF and >80% in ISF).

BMNC Culture in SF and ISF
BMNC were plated in 24 well culture plates (2x10 6 viable cells/50 µL DPBS/well) and covered with 500 µL SF or ISF. Cell viability was assessed at baseline using trypan blue and ranged from 74-96% across horses. Well contents were carefully mixed and plates incubated at 37°C in 5% CO 2 and 90% humidity. Remaining SF and ISF was preserved at 4°C for later addition of medium to replenish cell nutrients (200 µL added every 48 hours). All conditions and time points were performed in duplicate wells with cells from one well used for RNA-sequencing and the other for flow cytometry (macrophage activation markers CD14, CD86, CD206 and IL-10). Conditioned medium was aspirated at the same time points of cell harvest (48 and 96 hours and 6 and 10 days) centrifuged, and the cell-free supernatant used for cytokine and growth factor quantification (FGF-2, GM-CSF, IL-1b, IL-6, MCP-1, IL-10, TNF-a, SDF-1, IGF-1, IL-1ra, and PGE 2 ) using a PGE 2 ELISA kit (KGE004B; R&D Systems) and the Milliplex Map Equine chemokine/cytokine bead based array (Eqcttmag-93K,; MilliporeSigma). Details and findings from flow cytometry and cytokine and growth factor quantification are reported elsewhere (50) and were used in this study as a trait for weighted gene co-relation network analysis (WGCNA).

RNA Isolation and Sequencing
Cultured cells were recovered in 10 mM EDTA, centrifuged (12,

Functional Genomics
We adopted a multidisciplinary approach to functional genomics by employing several bioinformatics tools to tease out the biological significance of our data. We used WGCNA and DAVID in a semi-supervised analysis to identify biological processes of interest, and IPA to identify upstream regulators and activated and inhibited signaling pathways. By using this approach, we took advantage of both the superior annotation of biological processes from DAVID and the better pathway annotation of IPA. Together, these tools enabled us to make associations to our previous clinical studies to start to draw clinical translations to our findings. Weighted gene co-relation analysis was performed using WGCNA version 1.66 package in R to construct gene coexpression networks as described elsewhere (58,59). Gene coexpression clusters were generated from the whole transcriptome in SF and ISF datasets separately over time. Only genes expressed in at least 50% of samples in each dataset were included in the analysis (16,318 genes in SF and 18,038 genes in ISF). In order to normalize the data, the GeTMM values for each gene were log2 transformed. Next, a pairwise correlation matrix was constructed between all pairs of genes across the samples, and a matrix of weighted adjacency was generated by raising co-expression to a power b = 9, as determined for our sample set (58,60). A topological overlap matrix (TOM) was then assembled and used as input for hierarchical clustering analysis. Then, a dynamic tree cutting algorithm was used to identify gene clusters or modules (i.e., genes with high topological overlap) in an unsupervised fashion. Gene modules were visualized by heatmap plot (TOMplot) of the gene network topological overlap. Module relationships were summarized by a hierarchical clustering dendrogram and TOMplot of module eigengenes (MEs). Associations between gene modules and traits of interest were tested by correlating MEs to trait score. Module-trait correlations were visualized using a heatmap plot and only modules with trait relationship significance (R 2 ) higher than 0.7 and a p-value ≤0.05 were considered for further analysis. Traits of interest used for WGCNA included: timeline from our previous study, previously reported CD14, CD86, CD206 and IL-10 expression measured by flow cytometry, and IL-10, IGF-1, MCP-1, IL-1b, TNF-a, PGE 2 and SDF-1 concentrations quantified in conditioned SF and ISF (50). Module memberships (MM; correlation between each gene expression profile (GeTMM) and the ME of a given module as an indicator of the intramodular connectivity) and gene significance (GS; correlation between the gene expression profile (GeTMM) and the trait score (e.g. cytokine concentration in conditioned SF/ ISF) as a measure of biological relevance) were calculated (58). Genes (network nodes) having MM ≥ 0.90, P-value < 0.05, and GS ≥ 0.5 were identified as intramodular hub genes (61). Gene ontology (GO) analysis was performed on the entire gene list derived from each module as described above using DAVID Bioinformatics Resources version 6.8 (62) to functionally annotate their biological processes (BP). Of note, no single time point was chosen to determine the module-trait correlations. The entire timeline of the study was itself a trait. Therefore, genes within each module were co-expressed at all time points and thus dominant (overrepresented) BPs for a given module were the same at all time points.
To predict upstream regulators relevant for each set of DEGs, analysis was performed using the Ingenuity Pathway Analysis software (IPA, 2018) (63). The analysis output provided a P-value of overlap, activation Z-scores, and the downstream targets for each predicted upstream regulator. Z-scores were used to predict activation state (activation or inhibition) of each upstream regulator/signaling pathway. Predicted upstream regulators were considered significant if they had P < 0.05 and activation Z-score >2 (activated) or <−2 (inhibited). Subsequently, we investigated overlap between the predicted upstream regulators for each set and the DEGs from the same set to identify potential regulators among those DEGs. Genes in common between the two analyses with Z-scores (generated by IPA) matching the direction of fold change (generated by DESeq2) were defined as potential regulators. To investigate the interaction and relationships between potential upstream regulators, all known proteinprotein interactions were referenced and matched using STRING version 10.5 (64). Potential upstream regulators of high interaction were selected to have their protein expression assessed in synovium from inflamed joints treated with BMNC or DPBS, as a means of identifying candidate biomarkers of BMNC-mediated resolution. Synovial membrane samples were obtained from a related study in which BMNC therapy induced marked inflammation resolution (49) and represented inflamed joints treated with autologous BMNC or DPBS. IPA was also used to determine activated and inactivated signaling pathways, considering significance at P < 0.05 and activation Z-score >2 (activated) or <−2 (inhibited) and a -log(p-value > 1.3, which corresponds to p>0.05). For cases in which a large list of pathways met this criterion, those with a -log(p-value > 3 (FDR <0.01) were given priority attention.

Differential Gene Expression
Principal component analysis (PCA) of DEGs between BMNC cultured in SF and ISF showed clear differences in clustering patterns as early as 48 hours, and progressively diverged over time, representative of differences in BMNC response to normal and inflammatory environments ( Figure 2A). Volcano plots depict up and downregulated DEGs when comparing ISF to SF cultures at each time point ( Figure 2B). Vertical and horizontal comparisons were made, with the number of DEGs between any two conditions reported ( Figure 2C). In vertical comparisons, ISF cultures were compared to their SF counterparts at each time point. Counts of upregulated genes were most remarkable at 96 hours and 6 days. In horizontal comparisons, each subset of BMNC was compared with its nearest time point to analyze gene expression variation from BMNC response along the timeline.
From baseline to 48 hours, when myeloid progenitors in BMNC commit to the monocyte/macrophage lineage, the number of DEGs was highest among all time points, for both SF and ISF. The expression patterns of DEGs identified in all vertical and horizontal comparisons were also visualized by heat map ( Figure 3A). Further, Venn diagrams were used to illustrate the intersection between DEGs identified by horizontal comparisons and revealed that SF and ISF cultures shared 64.9% of DEGs over the 10 days, while those expressed exclusively in the SF or ISF dataset represented 15.0% and 20.1%, respectively ( Figure 3B). Upset plots elucidating the intersection between DEGs identified by vertical comparisons revealed that most DEGs were exclusively expressed in ISF cultures at 6 days ( Figure 3C). An entire list of DEGs is available at Supplementary Table 1. interaction or connectivity (hub genes) among SF and ISF datasets separately (hub genes are denoted by bold cells in Supplementary Table 2). Co-expression analysis of 18,038 genes in ISF identified 11 module eigengenes (i.e., clusters) ( Figure 4A). Among these, modules turquoise, green, blue, brown, black and pink were positively associated with three of the assigned traits. The turquoise and green modules were positively associated with IL-1b concentrations in conditioned ISF, and thus interpreted as having an overall pro-inflammatory nature. Modules blue, brown, black and pink were positively associated with the timeline. The blue module was also associated with CD86 expression assessed by flow cytometry, which denotes macrophage activation (50). In the SF dataset, analysis of the 16,318 genes identified the same 11 gene modules; however, the timeline was the only trait with a positive relationship to MEs, which like ISF included the blue, brown, black and pink modules ( Figure 4B). Thus, the IL-1b-related green and turquoise modules are the standout, inflammation-associated differences between ISF and SF cultures, while the blue module also differed in the number of positively associated traits. Since events associated with inflammation resolution would only be present in an inflammatory environment, further dissection of WGCNA findings were centered on data from ISF cultures, while data from SF cultures were used as a point of comparison.

Co-Expression Network Analysis From BMNC in Response to Inflammation
To assess the temporal behavior of each module, the mean expression profile (mean GeTMM values for all genes) for each module was plotted over time for both SF and ISF separately ( Figure 5A). This comparison revealed that the blue and brown modules exhibited increasing mean expression profiles that similarly dominated over time in both SF and ISF, which however, were higher in ISF. The pink and black modules had lesser expression among modules identified which completely overlapped between SF and ISF, and therefore are not graphically represented. Additional comparisons for a given module between ISF and its SF counterpart also included the functional annotation of the genes within such modules ( Figures 5B-E).
The blue, brown, black and pink modules completely overlapped between ISF and SF regarding their gene list, HUB genes list (Supplementary Table 2) and functional annotations (Supplementary Tables 3, 4). Exclusive to ISF, the green module peaked at 48 hours, the same time at which the blue and brown modules started to exhibit increased expression in comparison to its SF counterpart. The mean expression profile in the ISF's turquoise module progressively decreased from baseline. Functional annotation of genes within each module was then inspected (Figures 5B-E).

Modular Gene Ontology Enrichment and Overrepresented Biological Processes (BPs)
Overrepresented BPs were ranked based on fold enrichment and having an FDR <0.05 (Supplementary Table 3). For cases in which a large list of BPs met this criterion, BPs with FDR <0.01 were given priority attention. Fold enrichment was determined by comparing the background frequency of total genes annotated to a certain BP in the specified species to the sample frequency of genes under such BP. Overrepresentation was defined by a positive fold enrichment value (65). Since the gene list for modules blue, brown, pink and black completely overlapped between SF and ISF (Supplementary Table 2), overrepresented BPs in any of these modules were the same for both groups ( Figures 5B-E). For the blue module, isoprenoid biosynthesis was the most overrepresented of the 79 BPs identified by GO.
Given the high number of overall (n=4148) and hub genes (n=699) in this module, a diversity of BPs was identified within it, and is collectively discussed below. Of note, genes identified with pro-resolving functions in the related previous studies (IL10 and IGF1) (49, 50) also allocated to the blue module. In the brown module, while the most overrepresented BP was "antigen presentation via MHC class II", most BPs in this module related to mitochondrial response to oxidative stress and energy metabolism homeostasis. Overrepresented BPs in the pink and black modules constituted a minor list and were associated with a variety of cell homeostasis and housekeeping functions. Exclusive to ISF cultures, BPs in the IL-1b-associated green module were primarily associated with macrophage response to damage, including mitosis, adjustment of lipid and glucose metabolism following circadian distress, activation of the amphireguling-STAT3 axis (GO:0032355~response to estradiol) and noncanonical NF-kB signaling, thus, likely a module with a pro-inflammatory signature. In the turquoise module, also associated with IL-1b production, overrepresented BPs reflected the response of myeloid progenitors to stress and IL-4 signaling, a key event in the response of pro-resolving macrophages, amplifying chromatin opening for enhanced mRNA transcription (44)(45)(46). In summary, ISF triggered an early pro-inflammatory response in BMNC progenitors (green module) leading to macrophage commitment and priming (turquoise module). These events enhanced the constitutive expression of homeostatic mechanisms from macrophages (blue and brown modules) required to counteract damage and recover homeostasis (Figures 5A-E).

Pathway Analysis, Upstream Regulators and Their Network Interactions
Ingenuity Pathway Analysis revealed activated and inactivated pathways in SF and ISF cultures ( Table 1, Supplementary  Table 4). Our pathway analysis results from 0-48 hours (performed with IPA) agrees with findings from GO analysis and points repeatedly to activation of the mevalonate pathway and isoprenoid biosynthesis (superpathways of cholesterol biosynthesis, geranylgeranyl diphosphate biosynthesis, cholesterol biosynthesis I, II and III, and mevalonate pathway I). The patterns of expression of genes involved in these pathways, comparing BMNC cultured in SF and ISF (Figure 6), highlight the increased expression of genes such as ACAA2, HADHA ACAT2 and FDPS in ISF, essential for the synthesis of isoprenoids and mitochondria beta-oxidation of fatty acids. Additional pathways activated at 0-48 hours included unfolded protein response in agreement with overrepresented BPs in the blue module, and estrogen biosynthesis, in agreement with the BP "response to estradiol" from the green module, peaking at 48 hours and progressively decreasing.
From 48-96 hours, the PPAR-signaling pathways were repeatedly identified, which agrees with the identification of PPARG as highly connected upstream regulator, depicting the activation of the PPAR-g and PPAR-a signaling pathways. Comparisons for BMNC cultured in SF and ISF for the expression of genes involved in these pathways ( Figure 6), highlight the higher expression of PPAR genes in ISF. It also highlights the higher expression of NFKB2 and RELB genes, encoding for drivers of non-canonical NF-kB signaling, which has essential pro-resolving functions. After 96 hours of culture in ISF, a mix of pathways involved in cartilage metabolism (Heparan Sulfate Biosynthesis), cellular homeostasis/ inflammation resolution (Cell Cycle: G2/M DNA Damage Checkpoint Regulation, Unfolded Protein Response) and leukocyte migration during inflammation (Netrin Signaling) were observed to be activated. This mixed profile may have resulted from continuously challenging BMNC with ISF every 48 hours as performed in our model (50). Inhibited pathways included leukocyte extravasation signaling, IL-6 across the time course, IL-15 production and acute phase response signaling, key players in the development and maintenance of synovitis and degenerative processes observed in osteoarthritis. While the "Osteoarthritis Pathway" was identified as the third most activated pathway in SF at 6-10 days, neither the model used in our study, nor the list of genes involved in such functional annotation support such a finding.
Twenty-three potential upstream regulators were identified as activated (p <0.05 and a Z-score ≥2) in ISF and 35 in SF ( Figure 7A, Supplementary Table 5, Supplementary  Figure 1). Within these, 5 transcription factors (EIF4E, LARP1, MAFB, NFE2L2 and SIRT2 genes), the transmembrane receptor TREM2, the enzymes LPL and PIK3R1, and the multifunctional receptor GABARAP were conserved between ISF and SF. Inactivated upstream regulators (p <0.05 and a Z-score ≤ 2) were also identified in both ISF (n=17) and SF (n=32). Within these, IL1B and its downstream signaling transcription factor RELA, the mitochondria fission receptor MFN2, and the transcription factor GATA1 were conserved between ISF and SF.

Tissue Expression of Upstream Regulators and Genes Central to Activated Pathways
From the patterns of gene expression for each activated upstream regulator assessed over time (GeTMM counts and fold change, Supplementary Table 6), those with a high network interaction and exhibiting evident fold changes over time (decreasing or increasing) were selected for immunohistochemistry. Genes identified as essential drivers of the most activated pathways, were also selected as IHC targets. Ultimately, available immunohistochemistry targets with antibodies known to crossreact with equine proteins included PPARg, phospho-PPARg, PPARGC1A, MVK, HMGCS1, IL-1b, MAFB, SIRT2, and CSF1. Tissue expression for phosphorylated PPARg (p=0.0032) and mevalonate kinase (p=0.0291) were lower in BMNC-treated joints (Figure 8), likely because the resolution process had already been achieved, and is similar to the patterns of expression of IL-10 in our previous study using the same synovial samples (49). None of the remaining histochemical findings reached statistical significance. Tissue expression for PPARGC1A was consistently higher in BMNC-treated joints for 5/6 horses compared to those treated with DPBS. In contrast, CSF1 and MAFB tended toward higher expression in DPBStreated joints. No differences were observed in the staining patterns for SIRT2 and IL-1b.

DISCUSSION
In this study, we identified differential transcriptional signatures of BMNC in response to ISF and SF. These same BMNC were previously shown to resolve synovitis following exposure to an inflamed synovial environment in vivo and in vitro (49,50). We demonstrate a temporal behavior of co-expressed gene networks and their association with traits from our previous studies (49,

50)
, as well as with the expression of key proteins in the synovium by immunohistochemistry. Our findings illustrate the elaborate balance of pro-and anti-inflammatory mechanisms shifting dominance through the recovery of joint homeostasis. BMNC responded to ISF with an early proinflammatory response (green module), characterized by a short spike in the expression of NF-ƙB-related genes, coincident with the peak of IL-1b secretion in conditioned ISF (50). This response was associated with increased expression of the blue and brown modules, 2 gene networks with homeostatic functions comprising known drivers of resolution, which were more highly expressed in ISF and were dominant among activated upstream regulators. Significant differences in the expression of phosphorylated PPARg and mevalonate kinase in synovial membranes from inflamed joints treated with BMNC, and equivalent IL-1b staining between BMNC-and DPBS-treated joints, emphasize the fine tuning of so-called pro-inflammatory pathways that must remain active at physiological levels during the resolution process. These observations highlight the differences between the pro-resolving effects associated with BMNC therapy compared to the anti-inflammatory effects observed following clinical treatment with corticosteroids (30, 49-51, 69, 70). The short spike of expression of the pro-inflammatory green module observed in ISF, is essential to trigger a cascade of events that culminates in a pro-resolving response. NF-ƙBrelated and mitogen genes co-expressed in the green module (Supplementary Table 3) play a key role in promoting the proliferation of macrophages necessary to counteract damage (24,25). NF-ƙB-related genes also induce increased expression of A B FIGURE 8 | Immunohistochemistry of synovial membranes from 6 horses with experimental synovitis treated with BMNC or DPBS. Marked resolution of inflammation was evident following BMNC therapy (31). Selected targets were potential upstream regulators of high network interaction that were differentially expressed in ISF cultures over consecutive time points, or central drivers of most activated pathways. (A) Scatterplots of composite staining scores (median, 95% confidence interval) for PPARg (PPARG gene), phosphorylated PPARg, PPARg co-activator 1 alpha (PPARGC1A gene), colony stimulating factor 1 (CSF1 gene), mevalonate kinase (MVK gene), 3-Hydroxy-3-Methylglutaryl-Coenzyme A Synthase 1 (HMGCS1 gene), interleukin-1b (IL-1b gene), sirtuin 2 (SIRT2 gene), and transcription factor MAFB (MAFB gene). Each dot in the scatterplot represents the composite score for each individual horse. Tissue expression for phosphorylated PPARg (p=0.0032) and mevalonate kinase (p=0.0291) were lower in BMNC-treated joints. (B) Representative sections of synovium from inflamed joints treated with BMNC or DPBS and stained for selected markers detailed above (scale bars 100µm).
genes with anabolic and anti-inflammatory functions within the blue and brown modules, which have critical roles in driving resolution of joint inflammation. The top four overrepresented BPs in the green module relate to a group of genes sharing important mitogenic activity. Genes such as MAD1L1 and CHAMP1 encode proteins that interact and regulate cell structure organization preceding mitosis (71,72). NR4A3 and NR1D2 encode transcriptional activators involved in proliferation, survival and differentiation of myeloid progenitor cells, and in adjusting myeloid progenitor cell metabolism to oxidative stress (73,74). NR1D2 does so by activation of IL-6 transcription, which is also required to induce expression of the IL-4 receptor and related downstream regulatory functions of macrophages, including their self-renewal (75)(76)(77). The BP "response to estradiol" was characterized by the expression of the amphiregulin (AREG) and STAT3 genes. Macrophages are an important source of amphiregulin produced during acute inflammation. The AREG/ERK/STAT3 signaling axis is required for the differentiation of progenitor cells during tissue repair and establishing a pro-resolving response (78)(79)(80)(81). AREG was more highly expressed in ISF and exhibited progressively decreasing expression over time, as the resolving response progressed (Supplementary Table 6). Additionally, estrogen accelerates the resolution of inflammation through the regulation of IL-10/STAT3-mediated deactivation of proinflammatory responses, such that post-menopausal women are prone to developing chronic inflammation (82). STAT3 was the activated upstream regulator in ISF cultures with the highest connectivity. The NIK/NF-ƙB signaling in the green module was highlighted by the expression of RELB and NFKB2. Both RELB and NFKB2 are subunits of the noncanonical NF-kB signaling, which in macrophages, can exert both pro-and anti-inflammatory effects (83,84). Non-canonical NF-kB signaling is critical to produce SDF-1a and recruit monocytes to the site of damage immediately following injury (83). Further, during monocyte-macrophage differentiation, non-canonical NF-kB signaling prevents hyperactivation of new macrophages by accelerating the removal of RelA and c-Rel (canonical NF-kB subunits) from pro-inflammatory gene promotors preventing overt inflammation. As such, blocking non-canonical NF-kB by inactivation of its IKKa subunit results in increased inflammation (84). In both SF and ISF cultures, RELB expression was positively regulated, with decreasing expression over time, while RELA was downregulated ( Figure 6; Supplementary Table 6). Combined, these signatures illustrate a fraction of molecular drivers of the acute response of BMNC to inflammation, which also sets the stage for establishing a pro-resolving response.
Increased expression of the blue module in response to inflammation, in parallel with the surge of the green module, and over a timeline associated with resolution in our previous studies (49,50), suggests a pro-resolving identity. Genes encoding established drivers of joint inflammation resolution (IL-10, IGF-1) allocated to the blue module (Supplementary Table 2). A major functional signature of this module was the activation of the mevalonate pathway and isoprenoids biosynthesis, comprised by the expression of genes encoding central drivers of the mevalonate/isoprenoid pathway (COQ2, HMGCR, FDPS, HMGCS1, GGPS1, MVK, PDSS1, PDSS2, GGDPS1, FDPS, ACAA2, HADHA and ACAT2. In agreement, the super pathway of cholesterol biosynthesis, geranylgeranyl biosynthesis (an isoprenoid) and mevalonate pathway had the highest activation scores in ISF cultures between 0 and 48 hours, as detected by IPA. The roles of the mevalonate pathway in steroidogenesis, counteracting oxidative stress and inflammation resolution, are well documented in the macrophage response to damage and inflammation resolution (85)(86)(87)(88). Deficiency of mevalonate kinase (MVK), a key enzyme in the mevalonate pathway, causes reduced synthesis of isoprenoids, leading to mitochondrial damage, subsequent oxidative stress and severe inflammation (89,90). Importantly, exogenous isoprenoid treatment in models of inflammation induces decreased oxidative stress and production of inflammatory markers, by increasing expression of the NF-kB inhibitor IkBa and antioxidant selenoproteins (86,(91)(92)(93). Gene expression for MVK and HMGCS1, key enzymes in the mevalonate pathway, exhibited a similar expression profile between themselves ( Figure 6; Supplementary Table 6). These enzyme genes were more highly expressed in SF, suggesting that inflammation negatively affects this pathway. A similar pattern of gene expression for IL10 and IGF1 was observed in our previous study, in which both genes were more highly expressed in SF than ISF (50). However, concentrations of IL-10 and IGF-1 in ISF were higher than in SF, because the decreased relative production in ISF was compensated for by higher macrophage counts (50). Our ongoing lipidomic study on BMNCconditioned SF and ISF from the same samples used in this study may help elucidate if a similar context applies for isoprenoid production. Interestingly, it was recently evidenced that transcriptional dysregulation of the mevalonate pathway is a key signature of the overt inflammation caused by SARS-CoV-2 infection, highlighting its importance for inflammation resolution (88). Differences in MVK expression detected in synovial membranes suggests that the in situ activity of MVK in synovitis resolution happens earlier in time, as suggested by our pathway analysis and in vivo study (49).
The second overrepresented BP in the blue module "proteasome ubiquitin-dependent protein catabolism" was comprised of genes that characterize the formation of the 26S proteasome. Inflammation-derived oxidative stress damages nascent proteins that become misfolded and targeted for degradation (94). The 26S proteasome is essential for the degradation of these proteins, preventing aggregate formation, which is part of the pathogenesis of several conditions (94). Such observation from GO analysis agrees with the IPA findings where the "Unfolded Protein Response" was identified among the most activated pathways following isoprenoid biosynthesis through the mevalonate pathway. The next overrepresented BP was "protein localization to the Cajal body". These are coiled bodies found in the nucleus of proliferating or metabolically active cells and are implicated in telomere homeostasis (95,96). This BP was represented primarily by genes encoding chaperonin containing tailless (CCT) proteins that are critical regulators of telomerase folding and trafficking (97). Depletion of CCT proteins cause Cajal body and telomerase mislocalization and failure of telomere elongation (97). CCTs are also required for folding of cytoskeletal proteins during cell proliferation (98). The fourth dominant BP, "response to thyroid hormone stimulus", reflects the effects of triiodothyronine on regulation of macrophage maturation and responses. Such responses include controlling cell migration and conferring protection against endotoxemia and LPS exposure, in great part through proliferation of tissue resident macrophages (99,100). In addition, half of the genes characterizing this BP encode for cathepsins that mediate endolysosomal protein degradation. In summary, the dissection of a minute part of the blue module, illustrates the effect of isoprenoids and thyroid hormones in improving metabolism and performance of BMNC-derived macrophages to counteract the effects of inflammatory oxidative stress (85,89,90,99,100). This response is, at least partially, achieved by adjusting proteostasis through the 26S proteasome and Cajal bodies, preventing degenerative protein aggregate formation during increased cell transcription, proliferation and metabolism in response to inflammation (94,97,98).
The signature of the brown module were BPs that comprised of a series of gene groups encoding proteins that regulate the mitochondrial respiratory chain and mitochondria-mediated regulation of energy metabolism (101,102). There is growing evidence of the pivotal role of mitochondria in energy metabolism adjustments required for inflammation resolution as shown in the brown module (103)(104)(105). In the face of inflammatory challenges, enhanced cell respiration induces oxidative stress, activating alternative sources of energy driving gluconeogenesis and enhancing mitochondrial fatty acid oxidation (101,102) as denoted by the increased expression of genes such ACAA2 and HADHA in ISF cultures ( Figure 6). These genes were, however, associated with the isoprenoid biosynthetic pathway, highlighting the recently reported role of mitochondrial isoprenoid biosynthetic process (106). Enhancing the mitochondrial respiratory chain is a homeostatic mechanism that prevents mitochondrial DNA damage and its subsequent cytosolic and extracellular release signaling through damage-associated molecular pattern (DAMP) receptors (107). The peroxisome proliferator-activated receptor-gamma (PPAR-g) co-activator 1a (PPRGC1A; PPRGC1A gene), a master regulator of mitochondria biogenesis and liver gluconeogenesis was among the outstanding genes involved in mitochondrial regulatory functions listed in the brown module (101) and exhibited a trend for higher expression in BMNC treated joints (Figure 8).
Following interaction with PPRGC1A, PPAR-g exhibits increased activity, interacting with a multitude of transcription factors and PPAR-g responsive elements (43). PPAR-signaling was an important activated pathway identified by IPA in both SF and ISF cultures between 48 and 96 hours (Table 1, Figure 6), denoting its homeostatic functions. Among different PPARs, PPAR-g was one of the most highly connected upstream regulators. In macrophages and other cells, PPAR-g signaling is a cornerstone of tissue repair and inflammation resolution, exhibiting myriad functions that are either PPAR-g-mediated or -dependent (43). Examples include shifting the production of pro-inflammatory cytokines towards anti-inflammatory and pro-resolving mediators, driving apoptosis and clearance of neutrophils, enhancing macrophage traffic, recruitment, phagocytosis and efferocytosis, improving mitochondrial respiratory performance, and the overall transcriptome of a regulatory response driving recovery of homeostasis (43,47,101,108). Of note, some isoprenoids, the signature of the blue module, as well as some specialized pro-resolving molecules, signal through PPAR-g, conferring increased production of IL-10, resistance to inflammatory stimuli and attenuated NF-ƙB activation following LPS stimulation (43,81,(109)(110)(111). Gene expression for PPARG was higher in ISF compared to SF. Significantly lower staining for phospho-PPAR-g expression in synovial membranes detected by IHC, combined with the timing at which PPARG was identified as an Upstream Regulator, suggest that its activity modeling synovitis resolution almost overlap with the acute phase of inflammation. Importantly, PPAR-signaling findings from IPA at 48-96 coincide with the time at which pro-resolving effects were observed in our previous in vivo and in vitro studies (49,50). Histochemical findings for phospho-PPAR-g and PPRGC1A in BMNC compared to DPBStreated controls suggests PPAR-g signaling was not only a BMNC response to the inflamed synovial environment, but also part of the beneficial effects of BMNC on treated joints. Moreover, there is recent evidence that PPRGC1A expression is required for chondrocyte metabolism and cartilage homeostasis, with PPRGC1A knockouts exhibiting delayed endochondral ossification, disruption of physeal morphology and severe premature osteoarthritis (112).
Also in the brown module, CSF1 was identified as an upstream regulator gene and was more highly expressed in ISF than SF cultures. The marked proliferation of BMNC in ISF observed in vitro (50) may also be a response from CSF-1 signaling (50). Following a spike of proliferation, stimulation of macrophage progenitors with interferon-gamma (IFN-g) or LPS induce maturation and cell cycle arrest, increasing MHC-II expression and developing the capability to quickly respond to inflammatory stimuli and antigen presentation (113). Given that synovitis in our model was induced by injection of LPS, it is not surprising that "antigen processing and presentation of polysaccharide antigen via MHC class II" was a dominant BP in the brown module, particularly considering that IFN-g was not detected by immunoassay in conditioned SF or ISF in a preliminary screening performed in our study (50). Combined, these observations from the brown module highlight the importance of macrophage proliferation and maturation, and PPRGC1A/PPAR-g signaling during macrophage-mediated joint homeostasis, identifying the need for further investigation of the therapeutic roles of PPAR-g-agonists in the recovery of joint health.
The combined overrepresented BPs in the turquoise module, exclusive to ISF, reflect the dynamic of serine/threonine phosphorylation, autophosphorylation and dephosphorylation during cytokine signaling and subsequent mRNA transduction (114). In our study, these events happened in response to IL-4 (50). IL-4 signaling regulates the establishment of a pro-resolving response in macrophages by adjusting chromatin conformation and access for RNA transcription initiation (45). It also represses the expression of classical pro-inflammatory genes, decreasing inflammasome activation, IL-1b production and pyroptosis (44). Alternatively, IL-4 signaling enhances DNA binding of RNAPII-pS2 and RNAPII-pS5, and H3K27Ac as a positive epigenetic modification that leads to increased expression of several gene networks with a pro-resolving function (44). STAT6, which also allocates to the turquoise module, is a major regulator of IL-4 (likely through IL-4r expression) and PPAR-g signaling (46). STAT6 facilitates binding of PPAR-g to DNA, which, like IL-4, increases the number of regulated genes and the magnitude of response from macrophages (44,46,47,115). IL-4r expression was increased in BMNC cultured in both SF and ISF, but was significantly higher in ISF (50), in agreement with the overall findings in this module. We did not, however, detect gene expression for IL-4 or the IL-4r agonist IL-13. Also, neither IL-4 nor IL-13 were detected in the synovial fluid (SF or ISF) in our in vitro and in vivo studies (29,31). Collectively, these findings suggest that, if produced by cells in the joint other than BMNC, their concentrations were either below detectable limits or absent, raising the possibility of signaling by an unidentified IL4r agonist or an unconsidered pathway with equivalent gene ontology. The decreasing mean expression of the turquoise module in relationship to the increasing mean expression of the blue and brown modules, indicates that the biological processes included in the turquoise module are primarily required for priming the macrophages in BMNC towards a pro-resolving response (45).
IL1B gene expression in BMNC and protein expression in synovial fluid coincided in time and duration, peaking at 24 hours and progressively decreasing over time (49,50). The IL1B gene was identified as a centerpiece among inhibited upstream regulators in ISF. Interestingly, IL-1b expression in the synovial membrane at 6 days following LPS model induction remained evident in both BMNC-and DPBS-treated joints. This finding suggests that the physiological levels of IL-1b expression required for homeostasis were conserved and not blocked, as commonly observed with the use of corticosteroids (51). Staining patterns for the other selected markers were inconclusive due to variability between individual horses. Although our experimental design was aimed at minimizing variability by including only horses of the same breed and a narrow range of age, variability is expected when working with human populations and animal models. Inbred mice poorly mimic the inflammatory reaction of people (116). While this study comprised a small cohort, the horse is an excellent model for the translation of inflammation (117) and an established model for the study of degenerative and inflammatory joint disease (1,118). In agreement with our in vivo (49) and in vitro (50) findings using the same horses, the combined effects of gene networks comprised by each module depict a pro-resolving response, dissecting some important signatures of this process, and identifying candidate biomarkers of synovitis resolution, such as PPAR-g and MVK expression and synovial fluid quantification of isoprenoids. The overall agreement between different analytical tools used in this study (between gene ontology of WGCNA-derived modules with IPA pathways and upstream regulator analysis) reinforce the significance of our data. However, future mechanistic studies are necessary to further determine the specific roles of pathways and biological processes identified in our study. A study of the transcriptome analysis from synovial samples from naturally inflamed joints treated with BMNC and DPBS would complement our current observations and further illustrate the molecular drivers of the synovial response to BMNC injection.

CONCLUSION
Our current data suggest that BMNC-derived mechanisms of resolution are primarily represented by constitutively expressed homeostatic mechanisms, whose expression is enhanced to counteract tissue damage. These homeostatic mechanisms translate into macrophage proliferation, enlarging the "macrophage army" to fight aggressors, also improving their general and mitochondrial metabolism to better resist the challenges of inflammatory oxidative stress. Such effect is partially achieved through the synthesis and signaling of lipid mediators that promote recovery of homeostasis. Further exploration of BPs and pathways not dissected in this study may identify additional targets for future investigations. The combined findings of our equine studies (27,49,50) and human clinical trials (48,119) highlight the long-lasting and superior pro-resolving effects of BMNC in the treatment of arthritic conditions. This study reveals important transcriptional signatures of BMNCinduced resolution of synovitis and reinforce that pro-resolving macrophages do not fit within commonly described pro-or antiinflammatory phenotypes established in artificial in-vitro systems (120,121). Current knowledge, including our study, suggests that in vivo, macrophages are by default homeostatic cells that, following injury, drive inflammation with the purpose of counteracting tissue aggressors, further guiding inflammation resolution and recovery of homeostasis (24,25,39,(122)(123)(124)(125). Our study also highlights candidate mechanisms by which BMNC provide lasting improvement in patients with OA. Therapeutic enhancement of PPAR-g signaling in joints with chronic inflammation may represent a novel strategy for resolving joint inflammation. Defining multiple mechanisms of macrophagemediated synovitis resolution may provide means to develop pharmacological pro-resolving therapies, bypassing the need for more invasive bone marrow aspiration and further advancing the treatment of many inflammatory arthropathies, not just OA.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author. The current RNAseq data was deposited in the Gene Expression Omnibus (GEO; GSE18552) repository.

ETHICS STATEMENT
This study was conducted in compliance with the Animal Welfare Act and the approval of the Virginia Tech Institutional Animal Care and Use Committee.

AUTHOR CONTRIBUTIONS
BCM, JNM, and LAD designed studies. BCM obtained and processed the samples. TSK and SCL mapped and quantified the sequenced libraries. HE-SA and SCL performed DEGs analysis. BCM and HE-SA performed functional genomics, figures preparation, data interpretation and conceptualization. KES performed immunohistochemical assays, which were scored by BCM, HE-SA, and KES. BCM prepared the manuscript. All authors edited, reviewed, and approved the manuscript.