Investigating mechanical and inflammatory pathological mechanisms in osteoarthritis using MSC-derived osteocyte-like cells in 3D

Introduction Changes to bone physiology play a central role in the development of osteoarthritis with the mechanosensing osteocyte releasing factors that drive disease progression. This study developed a humanised in vitro model to detect osteocyte responses to either interleukin-6, a driver of degeneration and bone remodelling in animal and human joint injury, or mechanical loading, to mimic osteoarthritis stimuli in joints. Methods Human MSC cells (Y201) were differentiated in 3-dimensional type I collagen gels in osteogenic media and osteocyte phenotype assessed by RTqPCR and immunostaining. Gels were subjected to a single pathophysiological load or stimulated with interleukin-6 with unloaded or unstimulated cells as controls. RNA was extracted 1-hour post-load and assessed by RNAseq. Markers of pain, bone remodelling, and inflammation were quantified by RT-qPCR and ELISA. Results Y201 cells embedded within 3D collagen gels assumed dendritic morphology and expressed mature osteocytes markers. Mechanical loading of the osteocyte model regulated 7564 genes (Padj p<0.05, 3026 down, 4538 up). 93% of the osteocyte transcriptome signature was expressed in the model with 38% of these genes mechanically regulated. Mechanically loaded osteocytes regulated 26% of gene ontology pathways linked to OA pain, 40% reflecting bone remodelling and 27% representing inflammation. Load regulated genes associated with osteopetrosis, osteoporosis and osteoarthritis. 42% of effector genes in a genome-wide association study meta-analysis were mechanically regulated by osteocytes with 10 genes representing potential druggable targets. Interleukin-6 stimulation of osteocytes at concentrations reported in human synovial fluids from patients with OA or following knee injury, regulated similar readouts to mechanical loading including markers of pain, bone remodelling, and inflammation. Discussion We have developed a reproducible model of human osteocyte like cells that express >90% of the genes in the osteocyte transcriptome signature. Mechanical loading and inflammatory stimulation regulated genes and proteins implicated in osteoarthritis symptoms of pain as well as inflammation and degeneration underlying disease progression. Nearly half of the genes classified as ‘effectors’ in GWAS were mechanically regulated in this model. This model will be useful in identifying new mechanisms underlying bone and joint pathologies and testing drugs targeting those mechanisms.


Introduction
Abnormal joint loading through skeletal malalignment, age and obesity are key risk factors for osteoarthritis (OA) (1).At least 12% of the OA population are younger, with post-traumatic OA (PTOA) arising from prior joint injury (2), and although progression rate varies, injury can cause debilitating chronic pain and reduced mobility within 10 years.Changes to bone physiology play a central role in the development of OA (3) manifesting clinically as disruption of the tidemark, subchondral bone sclerosis and osteophyte formation alongside articular cartilage destruction.These changes are mediated by bone resorbing osteoclasts, bone forming osteoblasts, and osteocytes whose primary role is to maintain the integrity and function of bone (4) in response to its metabolic, mechanical and inflammatory environment, by regulating gene and protein expression.The osteocytes, which comprise 90-95% of all bone cells, release a variety of factors in response to mechanical load including nitric oxide (NO), prostaglandin E2 (PGE2), and receptor activator of nuclear factor kappa-B ligand (RANKL) to regulate osteoclast and osteoblast function as well as acting in an endocrine manner, releasing factors that target distant cells in other tissues (5)(6)(7)(8).The extensive osteocytic network embedded within a lacunocanalicular system throughout the bone, with cell-cell and cellmatrix connections that extend to the bone surface (9,10) is ideally suited to the role of osteocytes as principal regulators of bone mechanosensation and mechanotransduction [reviewed in (11)].
In knee OA, mechanical loading causes pain, but the mechanism underlying the link between load and pain is unknown; bone is highly innervated and is one of the only tissues in the joint where structural changes correlate to pain in OA (12).Glutamate, the major excitatory neurotransmitter in the nervous system, also signals in peripheral and non-neuronal tissue.Both the regulation of glutamate release and the ability to respond to glutamate by expression of glutamate receptors (GluRs) has been reported in a range of joint cells including osteocytes (13)(14)(15).Glutamate signalling in bone is regulated by mechanical load (16) and linked to joint pain in humans (17).Increases in synovial fluid glutamate concentrations occur in both rheumatoid arthritis (RA) and OA patients (18) and correlate with an increase in inflammatory mediators in RA (19).In addition, synovial fluid glutamate levels in patients following anterior cruciate ligament (ACL) rupture and meniscal damage are comparable to that of OA patient synovial fluid levels and decrease with time post injury (13).Glutamate induces knee inflammation (20) and contributes to arthritic pain and swelling in an inflammatory arthritis model (21).Glutamate receptor antagonists within the joint alleviate symptoms of OA including pain, inflammation, and bone and cartilage pathology (13,22).Inhibition of AMPA/KA GluRs with NBQX at the time of onset, improves pain related behaviour in rat inflammatory arthritis (antigen induced) and mouse PTOA model (ACL rupture) and reduces swelling and inflammation (13,22).This protection is partly explained by NBQX inhibiting AMPA/KA GluR release of interleukin-6 (IL-6), an early driver of inflammation in both models of arthritis (15,22,23).IL-6 acts as a mechanosensitive cytokine playing a key role in the biochemical control of bone remodelling (24,25).Mechanical loading of osteocytes increases NO and IL-6 release and increases IL-6, osteoprotegerin (OPG), RANKL, and tumour necrosis factor-alpha (TNF-a) gene expression (26).Chronic IL-6 overexpression increases bone remodelling causing a net loss of bone by activating RANKLinduced bone resorption (27)(28)(29).Both pro-inflammatory and anti-inflammatory cytokines have been implicated in the pathogenesis of OA (30) with IL-6 among the most prominently elevated cytokine involved in the OA inflammatory response (31,32).Increased circulating IL-6, as well as increased body mass index, predicts development of radiographic knee OA (33) and single nucleotide polymorphisms in the IL-6 gene are associated with radiographic hand OA (34).In addition, serum IL-6 levels are strongly associated with the incidence of age-related OA (35) and a predictive marker for the risk of OA progression (36).The role of IL-6 in OA pathophysiology has been studied in OA animal models showing that IL-6 is largely destructive (37) affecting both cartilage and subchondral bone (38).
Although mechanical load is the major driver of human OA, humanised in vitro osteocyte models investigating mechanics and its interactions with inflammation are very limited and have not been validated against in vivo models or human clinical data (39,40).We have used our humanised in vitro model of osteocytes to detect early responses to either IL-6, a driver of degeneration and bone remodelling in animal (22, 37) and human (33, 34) joint injury, or mechanical loading, to mimic early OA stimuli in joints.Our osteocyte model (adapted from Vazquez et al. (41)) is derived from human Y201 stem cells differentiated into osteocytes in 3D type I collagen gels and mechanical strain applied using our custom loading device (adapted from Vazquez et al. (41)).The aim of this current study was to assess the effect of pathological mechanical load on the osteocyte signature and determine the influence of IL-6 on readouts that have been reported in OA.This will help identify mechanical and inflammatory mechanisms that cause pain or alter bone tissue structure in vitro and provide new mechanistic insight into disease progression.

Materials and methods
Chemicals were from Sigma (Poole, UK) and tissue culture and molecular biology reagents from Thermo Fisher Scientific (Invitrogen, Paisley, UK) unless otherwise stated and were of analytical grade or above.sIL6r and IL6 were from Peprotech, UK.

Cell culture
Y201 hTERT-MSCs, gifted from Prof Paul Genever (University of York), were used as a model human MSC line (42).These cells can rapidly differentiate in 3D to osteocytes in a manner similar to bone marrow MSCs in spheroid cultures (43) and maintain mechanoresponsive behaviour (44).Cells were cultured in basal medium [Dulbecco's modified Eagle's medium (DMEM), high glucose, pyruvate, GlutaMAX ™ , supplemented with 5% fetal bovine serum (FBS), 100 U/mL penicillin and 100 mg/mL streptomycin] at 37°C in a humidified atmosphere of 5% CO 2 , 95% air.At 80-85% confluency cells were sub-cultured by treating with TrypLE ™ .Y201 cells were incorporated into type I collagen gels at a concentration previously determined (41).Briefly, lyophilised rat tail tendon type I collagen was dissolved in 7mM glacial acetic acid and mixed 4:1 with 10X MEM containing 11g/L sodium bicarbonate on ice and neutralized [1M tris (hydroxymethyl)aminomethane (Tris) base, pH 11.5] to give 2mg/mL type I collagen gels.Y201 cells (0.125 x 10 6 cells/gel or 0.05 x 10 6 cells/gel for RNAseq) diluted in aMEM (<10% of total gel volume) were added to the collagen on ice and 250mL distributed into 48-well plastic (phenotype studies and IL-6 studies) or silicone (loading studies (41)) plates (Supplementary Methods 1), for polymerization at 37°C.After 1 hour, 800mL basal medium supplemented with osteogenic differentiation factors (50µg/mL ascorbate-2-phosphate, 5mM b-glycerophosphate, 1nM dexamethasone) was added onto the surface of the gels and cells cultured at 37°C with media changes every 3-4 days for the indicated periods.

Assessment of cell viability
Cultures grown in plastic plates for 7 days were rinsed with phosphate buffered saline, pH 7.3 (PBS), incubated with 1µL hoescht (1mg/mL) and 4µL propidium iodide (100µg/mL) in serum free medium for 2 hours at 4°C and then for a further 2.5 hours at 37°C before washing overnight at 37°C in normal culture medium with gentle agitation.Cells were fixed in 1% (wt/ vol) paraformaldehyde for 30 minutes at 4°C, washed in PBS prior to overnight infiltration with 50% OCT compound (Tissue Tek) in PBS at 4°C.Gels were frozen in fresh OCT compound onto cryostat stubs using dry ice and cryosections cut at 20mm using a Bright OTF5000 cryostat and collected on polysine slides (VWR, Lutterworth, UK).Slides containing sections were mounted in VECTASHIELD ® Mounting Medium containing DAPI (1.5 mg/ mL) to counterstain DNA (Vector Laboratories, Peterborough, UK) and viewed using a light microscope (BX61, Olympus).

Assessment of cell phenotype 2.1.2.1 Cell shape and immunolocalization of osteocyte markers
Cultures grown in plastic plates for 7 days were rinsed with PBS, and gels fixed in 1% (wt/vol) paraformaldehyde, frozen and sectioned as outlined above.Sections were stained with phalloidin to assess cell morphology or processed for immunocytochemistry.Phalloidin-iFluor conjugate staining was carried out according to manufacturer's protocol (Abcam).Slides were removed from -20°C, allowed to equilibrate to room temperature and hydrated with 2mL of PBS per section.PBS was aspirated off, sections permeabilised with 0.1% Triton X-100 in PBS for 5 minutes and washed 3x with PBS prior to treatment with 100mL of 1X phalloidin conjugate (1mL per 1mL of PBS + 1% BSA) per section.Sections were left to stain for 60 minutes at room temperature whilst protected from light and washed 3x with an excess volume of PBS prior to mounting in VECTASHIELD ® Mounting Medium containing DAPI (1.5 mg/ mL).Stained sections were visualised using fluorescent imaging (BX61, Olympus).For immunocytochemistry, each step was performed at room temperature unless stated otherwise and between each incubation step, sections were washed 3x 5 min in 0.01 M PBS containing 0.001% Tween 20 (wash buffer).All antibodies were diluted in wash buffer.Cells were washed before blocking in 2% (v/v) normal goat serum (Dako UK, Ely, UK) for 1 hour.After overnight incubation at 4°C with a rabbit polyclonal primary antibody to Sclerostin (Abcam; ab75914) diluted 1:100, cells were washed before incubating for 1 hour with goat anti-rabbit Alexa 488 conjugated secondary antibody (4 mg/mL; Molecular probes, Invitrogen).Finally, after washing, cells were mounted in VECTASHIELD ® Mounting Medium containing DAPI (1.5 mg/ mL).Representative cells from multiple fields of view were imaged by confocal microscopy (Leica TCSSP2, Germany) using a 63x oil immersion objective lens and appropriate settings for AlexaFluor 488 (green) and DAPI (blue).Negative controls where the primary antibody was omitted or replaced with rabbit IgG were devoid of fluorescent signal (data not shown).

RT-PCR analysis of osteocyte markers
RNA was extracted from cells grown in gels from duplicate cultures at day 7 and analysed for the expression of osteocyte markers by RT-qPCR.Briefly, gels were placed into 600µl buffer RLT Plus (RNeasy Plus kit, Qiagen) and 6µl b-mercaptoethanol added prior to loading onto a QIAshredder spin column to disrupt the gels.Samples were centrifuged for 2 minutes at full speed and the lysate added to gDNA Eliminator spin columns to digest genomic DNA.Samples were processed according to manufacturer's instructions (Qiagen) and RNA eluted in 30µl RNAse/DNAse free water.RNA quality and concentration were assessed by TapeStation Analysis (Agilent Technologies, UK). cDNA was generated in a 20 µL reaction from 500 ng RNA using 50 ng random hexamers (0.5 mg/mL; Promega) RNasin ® RNase Inhibitor (40U), 5mM DTT, 0.5mM each dNTP, and Superscript IV reverse transcriptase (200 units).Gene expression was measured by SYBR green RT-PCR using the Platinum ™ SYBR ™ Green qPCR mix and the AriaMX qPCR system according to manufacturer's instructions (Agilent Technologies UK) with 200 nM forward and reverse primers (Supplementary Table 1A) and the following cycle conditions: 1 cycle of 95°C, 3 minutes; 40 cycles of 95°C, 15 seconds and 60°C, 30 seconds; 1 melt cycle of 95°C, 1 minute, 65°C 30 seconds, 95°C 1 minute.

Mechanical loading of collagen gels 2.2.1 Strain validation in collagen gels
Force strain relationships were validated in collagen gels within the custom-built Cardiff loading device which comprises a deformable silicone multiwell plate within a 3D printed loading device (Supplementary Figures 1A, B).Defined vertical displacements applied with a Bose Electroforce 3200 machine (TE Instruments) extend the levers outwards and caused the plate to stretch.Collagen gels (2mg/mL) containing 500mL of blue-and violet-coloured microspheres (10mm, Polysciences, Park Scientific Ltd, UK) subjected to vertical displacements of 0-2.1mm at 0.35mm intervals were imaged by light microscopy and a tracking code, written in Matlab calculated the strain/displacements relationships (Supplementary Figures 1C-E).There was a linear relationship between displacement and strain up to displacements of 0.7mm.A displacement of 0.7mm represented a pathophysiological load of 4300µe ( ± 103) and was used for all experiments.

Mechanical loading
For loading, 3D Y201 cultures were prepared and cultured in the silicone plate in 800mL of basal medium containing osteogenic differentiation factors and incubated at 37°C in 5% CO2/95% air atmosphere for 5-days.After this time, the media were replenished and left for 24 hours.One hour prior to loading, media were removed and 800mL osteogenic media added.An hour later, silicone plates were loaded using a BOSE ElectroForce ® 3200 loading instrument (TE Instruments, UK) to stretch the plate causing cyclic compression in all wells (pathophysiological load 4300me induced by 0.7mm displacement, 10Hz, 3000 cycles (41,45,46).The loading regime was chosen to recapitulate in vivo models where validated high physiological strains induced osteogenesis (45,47,48).This high (pathophysiological) strain down regulated sclerostin (Supplementary Figure 1F) a known osteocyte derived mechanoresponsive molecule that is a potent regulator of bone formation (46) as well as inflammatory mediators relevant to osteoarthritis (41).In addition, osteocytes respond to mechanical load within seconds, revealing gene expression changes within 1 hour (49) and protein changes 24-72 hours later (41).Control gels in the silicone plate were placed into the loading device but received no load.Loading was controlled using WinTest ® Software 4.1 with TuneIQ control optimization (BOSE).Media was collected after 1hr and 24hrs, aliquoted and frozen (-20°C) for analysis of released factors.Gels, 1hr post load, were placed into 600µL buffer RLT Plus (RNeasy Plus kit, Qiagen) and stored at -80°C prior to RNA extraction.

Confirmation of the osteocyte response to load
RNA was extracted 1 hour post load using RNeasy Plus kits, RNA eluted in 30µl RNAse/DNAse free water, cDNA synthesised and RTqPCR performed as described above (section 2.1.2.1).Reference genes, 36B4, YWHAZ, RPL13A, 18S, b-actin, GAPDH were tested across experimental conditions.The geometric mean of YWHAZ and 18S (stability value 0.292), were identified by RefFinder (50) as the most stable and used to calculate fold change relative to untreated cells using the DDCT method (51).All primers (Supplementary Table 1A) were purchased from MWG and validated using a standard curve of five serial cDNA dilutions with primer efficiencies between 90-110% (52).

RNAseq analysis: generation of the osteocyte 'mechanosome'
RNA was extracted from loaded samples 1 hour post load (n=6) and unloaded controls (n=5) using RNeasy Plus kits as described above (section 2.1.2.2).RNA was eluted in 30µl RNAse/DNAse free water and RNA quality, and concentration assessed by TapeStation Analysis (Agilent).An RNA sequencing library was prepared for the mechanically loaded and control samples, using the New England Biolabs Ultra II directional RNA library prep kit (Wales Gene Park).cDNA was synthesized using this RNA which, after undergoing fragmentation, had adaptors ligated to the ends.The MiSeq Nano system (Illumina) was used to complete a sequencing library quality control after which sequencing was performed using the NovaSeq 6000 system (Illumina) running a 2 x 100bp pairedend reads run on a NovaSeq S1 flow cell.Trimming to remove adapter sequencer and poor-quality ends of reads was performed by Trim Galore using default parameters in paired-end mode.Trimmed paired-end reads were aligned to the GRCh38 no_alt_plus_hs38d1 analysis set reference using STAR (v2.5.1b), an ultrafast universal RNAseq aligner, following the 2-pass method (53).QC metrics were generated using FastQC (v0.11.2), and summary statistics were generated using Samtools (v0.1.19)flagstat.Raw counts were calculated for all samples for both (i) exons and (ii) genes using Subread featureCounts Version 1.5.1.Counts were generated for paired end read fragments summarized at exon level and then aggregated at transcript level.Ambiguity between the labels on sample C1 and L1 and the suspicion that these 2 samples had been inadvertently switched led to an additional PCA analysis which showed C1 to cluster with loaded samples (Supplementary Figure 2); because of this C1 and L1 were excluded from further analysis leaving n=4 controls and n=5 loaded samples.Differentially expressed genes were identified using an DEseq2 analysis (54) on normalised count data.The resultant p-values were corrected for multiple testing and false discovery issues using the FDR method (55).RNAseq data are available from Mendeley Data (DOI: 10.17632/5md5rnybcs.1).

Human osteoprotegerin and glutamate ELISAs
Media collected 1 hour, and 24 hours post-load were analysed for the release of Glutamate (KA1909, Bio-Techne) and OPG (AB100617, Abcam) using commercial kits following manufacturer's instructions.

Cytokine analysis
Frozen media samples were thawed on ice and centrifuged for 5 minutes at 3000 rpm, and 23 cytokines profiled (Supplementary Table 1B) in each sample by Luminex bead based multiplex assay using a Merck Milliplex ® MAP human cytokine/chemokine magnetic bead panel kit (2923824 HCYTOMAG-60K-23) following the manufacturer's instructions.

sIL6r/IL6 treatment of Y201 cells in 3D gels
Y201 cells were embedded in 3D type I collagen gels as described above at a density of 0.125 x 10 6 cells/gel and grown in 48-well plastic plates for 24 hours in basal media at 37°C, 5% CO2.Media was replaced with 800mL osteogenic media and cells cultured for 7 days with media changes every 3 days.At day 7, the media was replaced with 800mL osteogenic media containing IL-6 (5ng/mL) and sIL-6r (40ng/mL) in all but the control wells (26, 31, 32).At 24 and 72 hours, media was removed and stored at -20°C for cytokine analysis.TRIzol ™ reagent (500mL) was added to each gel at 24 hours and pipetted repeatedly to dissolve the gels and lyse the cells, prior to storage at -80°C.

RNA extraction and RT-qPCR analysis of gene expression
Total RNA was extracted from gels using TRIzol ™ reagent according to the manufacturer's protocol.RNA was DNase treated to remove genomic DNA (Ambion; Applied Biosystems, UK) and re-suspended in 50 µL RNase-free water.RNA integrity and concentration were assessed by Nanodrop ™ .cDNA was generated in a 20 µL reaction from 150 ng RNA as described above (section 2.1.2.2).RT-qPCR was carried out as described above using the geometric mean of EEF and RPL13A (RefFinder stability value 0.375) for normalisation.

Cytokine analysis
Aliquoted media samples were centrifuged to remove cells and supernatants vortexed briefly prior to use.A multiplex electrochemiluminescence (ECL) kit (Supplementary Table 1C; U-Plex Proinflam Combo 1 Human K15049; Meso Scale Discovery, USA) and single-plex ELISAs (Glutamate KA 1909 Abnova, OPG RDR-OPG-Hu 2bScientific) were utilised to measure levels of released molecules according to manufacturer's instructions.Multiplex ECLs were carried out in the Central Biotechnological Services (Cardiff University) utilising a Mesoscale discovery (MSD) plate reader to determine chemiluminescence measurements.

Statistics and data analysis
Results are presented as mean ± SEM.Graphs show individual data points, box and whisker plots of minimum and maximum values, 25th and 75th quartiles and median.Data were tested for normality and equal variances prior to transformations where necessary and appropriate statistical testing as indicated in the figure legends (Minitab 20).Differences were considered significant at p=0.05.For all statistics, unless stated otherwise, treatments were compared to untreated controls.Functional gene enrichment and pathway enrichment analysis were performed on up and down DEG sets using gProfiler (56) and Enrichr databases (57-59), respectively.RNAseq data was compared to the osteocyte signature (60) as well as a list of gene ontology terms (GO terms), compiled using AmiGO (61)(62)(63) and the search terms 'bone', 'pain', 'inflammation', and 'mechanical load'.For phenotype studies, data is collected from 4 independent studies (N=4).For IL-6 effects on bone remodelling, OPG gene data is representative of 3 independent studies (N=3).

Differentially expressed genes, model reproducibility and the osteocyte signature
RNA extracted from cells differentiated in silicone plates 1 hour post load was of high quality (RIN scores >9; Supplementary Figures 3B).RT-qPCR analysis revealed SOST expression was decreased by loading thereby confirming that the model cells were responding appropriately to load (Supplementary Figure 1F).RNA was processed for RNAseq and the dataset generated analysed for the presence of genes involved in the development and maturation of osteocytes and the mineralisation process (Figure 2A).Our data was compared to publicly available data generated from IDGSW3 cells at day 3, 14 and 35 of differentiation (Figure 2B) (57), and osteocyte-isolated samples from the Osteocyte Enrichment cohort (Figure 2C) (65) and the relative temporal expression of various osteogenic markers during the transition from osteoblast to osteocyte summarized [Figure 2D; adapted from (10)].
Hierarchical clustering of the sample set grouped replicates and separated biological conditions (Supplementary Figure 2C).Principal component analysis (PCA) characterised further the experimental variability revealing loaded samples L2 and L3 to be different to L4-L6 (Supplementary Figures 2D, E).
In total 7564 genes were differentially regulated (Padj p<0.05): 3026 down and 4538 up regulated by mechanical load (Supplementary Table 2A).Of these, 3824 genes were up regulated and 532 down regulated fold change (FC) >2.A volcano plot representing the log of the adjusted P value as a function of the log ratio of differential expression shows differentially regulated genes as red dots or triangles, the latter corresponding to genes where log2 FC is too low/high to be displayed on the plot (Supplementary Figure 2F).We compared our in vitro osteocyte 'mechanosome' data to the 1004 protein encoding genes from the published osteocyte signature representing human genes enriched in osteocytes relative to bone marrow and other osteoblast lineage cells (60).This revealed that 937/1004 osteocyte signature genes were expressed (Supplementary Table 2B) and 67 were not expressed in our dataset (Supplementary Table 2C).Of these, 379 were regulated by mechanical load (248 UP, 131 DOWN; Padj<0.05;Supplementary Table 2D).Functional gene enrichment analysis of the DEGs using gProfiler revealed alterations in genes involved in several biological processes, molecular function, cell compartments, and the reactome (Supplementary Table 3A).Analysis of the up regulated gene set revealed enrichment of genes including those involved in metabolic processes, cell response to stress, and cell component organization as well as genes involved in post-translational modification and tolllike receptor signalling (Supplementary Table 3A GEA.UP).Down regulated genes were enriched in RNA processes, and cilium processes (Supplementary Table 3A GEA.DOWN).Additional gene ontology searches were performed for terms related to mechanical load including 'response to mechanical signalling', 'mechanically gated ion channels', 'mechanosensory behaviour', 'cell response to mechanical stimulus', and 'detection of mechanical stimulus' revealing several regulated genes (Supplementary Table 3B).Of the 166 genes listed across these terms, 46 genes were upregulated and 24 downregulated by mechanical load (padj <0.05).

Mechanical load of osteocytes regulates readouts of osteoarthritis
Gene ontology searches were performed for terms related to readouts of OA including pain (Supplementary Table 3C), bone remodelling (Supplementary Table 3D) and inflammation (Supplementary Table 3E).

Pain and the glutamate signalling pathway
Pathophysiological loading of osteocytes in our 3D model regulated genes known to be involved in gene ontology pathways linked to OA pain including genes related to sensory perception of pain and neuropathic pain and members of the glutamate signalling pathway (Supplementary Table 3C).Of the 253 genes listed across these terms, 43 genes were up regulated, and 22 genes downregulated by mechanical load (Padj < 0.05, Supplementary Table 2A).GRIK1, GRIA4, GRIN2B, GRIN2C, GRIN2D, GRIN3A, SLC1A2, SLC1A3, GRID2 were expressed but not regulated by load (data not shown).In addition, COL12A1 and COL16A1, present in bone marrow lesions (BMLs) (66) which correlate to OA pain (67), were upregulated 4-(padj=0.000)and 1.7-fold (padj= 0.049) by load, respectively (data not shown).

Bone remodelling
Pathophysiological loading of osteocytes in our 3D model regulated genes known to be involved in gene ontology pathways linked to bone remodelling including ossification, bone resorption, and bone mineralisation (Supplementary Table 3D).Of the 217 genes listed across these terms, 49 genes were upregulated and 31 downregulated by mechanical load (Padj < 0.05).In addition, RANKL (TNFSF11) and OPG (TNFRSF11B) were expressed but not regulated by mechanical load; SOST was not detected by RNAseq (data not shown).Pathway enrichment analysis of mechanically regulated genes using Enrichr (57-59) revealed associations with the bone remodelling/RANKL pathway  4).OPG and RANKL release was measured by ELISA.OPG release did not change over time in control cultures but was significantly increased by load at 24 hours (Figure 3B; 5.6-fold p=0.003; load at 24 hours vs load at 1 hour p=0.001).RANKL was below the level of detection in all cultures (data not shown).

Inflammation
Pathophysiological loading of osteocytes in our 3D model regulated genes known to be involved in pathways linked to inflammation including genes associated with an acute and chronic inflammatory response (Supplementary Table 3E).Of the 371 genes listed in these terms, 69 genes were upregulated, and 32 genes downregulated by mechanical load (Padj <0.05).Many of these genes also belonged to GO 3).Enrichr database analysis also revealed association with the NFkB pathway (Biocarta 2016: 12/ 21 p = 6.59E-04) (Supplementary Table 4).

IL-6 stimulation of osteocytes regulates readouts of osteoarthritis
Since IL-6 has been reported to be elevated in OA (31) and after knee injury significantly contributing to baseline KOOS and change in KOOS over 3 months (32), we stimulated our model with concentrations of IL-6 and its soluble receptor reported in synovial fluid from patients with OA (31) or following injury (32).

Pain and the glutamate signalling pathway
IL-6/sIL-6r treatment downregulated GRIA1 (3-fold, p=0.007) mRNA expression (Figure 5A).SLClA1 and SLC1A3 were expressed but levels did not change with IL-6 treatment (data not shown).Glutamate release into the media increased with time in culture in control (p=0.005) and IL-6/sIL-6r (p=0.006)treated cultures but there was no effect of IL-6/sIL-6r treatment on glutamate release (Supplementary Figure 4).

Bone remodelling
IL-6/sIL-6r treatment downregulated ALPL (2-fold; p=0.003),BGLAP (1.6-fold; p=0.042),TNFRSF11B (2-fold; p=0.029) mRNA expression (Figure 5B).Levels of Col1A1 did not change (data not shown).OPG release into the media increased with time in culture in control (p<0.001) and IL6/sIL6r (p<0.001) treated cultures but there was no effect of IL-6/sIL-6r treatment on OPG release (Supplementary Figure 4).The effect of mechanical load on cytokine release from osteocytes grown in 3D collagen gels.Media was analysed for the release of cytokines 1 and 24 hours post load using a Luminex bead based multiplex assay (Merck Milliplex ® MAP human cytokine/chemokine magnetic bead panel kit).Samples were compared to control at 1 hour unless stated otherwise (n=3/treatment; GLM ANOVA with Tukey post hoc tests; IL-6 and GM-CSF logged data; IP-10 and RANTES 2-sample t-tests).

A 3D model of human osteocyte like cells
Y201 cells embedded in 3D collagen gels for 7 days displayed appropriated dendritic morphology (68), and expressed the mediator of osteocyte mechano-responses, sclerostin (69).SOST is not expressed in the early stages of differentiation of the osteoblast lineage, but levels increase as the osteocyte matures and become surrounded by mineralized bone (70).In addition, several genes were expressed that have been identified as being involved in the development and maturation of osteocytes [reviewed in (71)(72)(73)].These included PDPN and CD44, markers of early osteocyte differentiation and required by osteocytes to initiate proper dendrite formation (74,75), and a number of cytoskeletal proteins involved in actin dynamics such as PLS3, which encodes for an actin bundling protein required for osteocyte cytoskeleton The effect of IL-6/sIL-6r on the expression of markers of (A) pain and (B) bone in osteocytes grown in 3D collagen gels.Gene expression was assessed at 24 hours post IL-6 (5ng/ml) and sIL-6r (40ng/ml) treatment and normalised to the geomean of two housekeeping genes (HKG), EEF and RPL13A.Samples were compared to control at 24 hours (One way ANOVA with Tukey post hoc tests; n=3).organization (76), destrin, CAPG, and CDC42 (77).GJA1 (CX43) which is required for osteocyte communication (78), and MMP14, which is involved in canaliculi formation (72,79) were also expressed along with DKK1, PHEX and FGF23 which play important roles in osteocyte development (72).BGLAP, COL1A1, TNFRSF11B, TNFSF11, PHOSPHO1, HDAC5, CYP19A1, RUNX2, DLX5, ATF4 and AP1 were also present representing genes that encode for proteins involved in osteocyte maturation and the mineralisation process [reviewed in (71)] (80).The genes involved in osteocyte development and maturation (5, 10, 70, 72, 74-79, 81-84) and found to be present in our dataset are summarized in Figure 7.Of the 1004 known human protein coding genes reported to reflect the in vivo osteocyte transcriptome signature (60), 93% were expressed in our model indicating that the cell phenotype and differentiation status is a good representation of osteocyte-like cells.The 7% of osteocyte signature genes not expressed in our RNAseq dataset included genes associated with the nervous system, the skeleton, angiogenesis, and cell function (Supplementary Table 2C).However, these may have been below the limit of detection since SOST was detected by RT-qPCR suggesting a reduced sensitivity in the RNAseq technology.

Mechanical loading regulates readouts of osteoarthritis in our 3D model of osteocytes
Mechanical loading decreased SOST expression in all 3D osteocyte cultures consistent with its regulation in vivo (46,85) and regulated 70 genes linked to known mechanical responses by GO term enrichment confirming the model's expected response to mechanical load.In total, mechanical loading of osteocytes down regulated 3026 genes and upregulated 4538 genes.Functional enrichment analysis using gProfiler revealed these were involved in many important cellular processes.Of note was the down regulation of processes involved in the primary cilia, an important key player in osteocyte mechanosensing (86) and bone formation in response to mechanical forces [reviewed in (87, 88)].Thirty eight percent of the previously reported osteocyte signature genes (60) were regulated by mechanical load in our 3D model consistent with the dogma that osteocytes are highly responsive to mechanical load (60).Load regulated several genes associated with skeletal diseases such as osteopetrosis, osteoporosis and osteoarthritis in keeping with the human orthologs associated with osteoarthritis and osteoporosis reported in the in vivo mouse The effect of sIL-6r/IL-6 treatment on cytokine release from osteocytes grown in 3D collagen gels.Media was analysed for the release of (A) IL12p70, (B) IL4, (C) TNF-alpha, and (D) IL2 24 and 72 hours post treatment with IL-6 (5ng/ml) and sIL-6r (40ng/ml) using a multiplex ECL kit (Meso Scale Discovery).Samples were compared to control at 24 hour unless stated otherwise (n=3-4/treatment; GLM ANOVA with Tukey post hoc tests).osteocyte transcriptome (60).Of the 25 genes of the osteocyte transcriptome signature associated with skeletal phenotype (60), 15 were expressed but not regulated in our model, 3 were upregulated (CADM1, KAZN, STARD13) and 3 downregulated (CC2D2A, LTBP1, PLS3) by mechanical load (Supplementary Table 5A).Interestingly, mutation in PLS3 which encodes for plastin-3 an actin bundling protein required for organisation of the cytoskeleton in osteocytes (76) results in X-linked osteoporosis (89).In addition, of the 211 genes enriched within the signature that are associated with a skeletal phenotype in the Mouse Genome Informatics database (MGI) (60), our model expressed 108 genes that were not regulated by load, 50 that were upregulated, and 37 that were downregulated (Supplementary Table 5B).Bone marrow lesions (BMLs) occur in areas of bone remodelling and are associated with subchondral bone microdamage and correlate with pain in OA (67,90,91).Their presence and location are associated with altered joint loading as occurs in joint malalignment and absence/regression of BMLs occurs following reduction of focal contact stress across the joint (92).Joint pain resolution is associated with diminished BMLs.Comparison of our mechanosome data to the 78 DEGS identified as BML hub (93) revealed 56 genes to be expressed in our dataset, with 8 down regulated and 24 upregulated by load (Supplementary Table 5C).
To determine whether mechanical load activated readouts classically associated with the structural and symptomatic changes in osteoarthritis, we performed gene ontology searches for pain, bone remodelling and inflammation and pathway enrichment analysis using Enrichr.

Pain
Mechanical loading of osteocytes in our 3D model regulated 26% of the 253 genes known to be involved in gene ontology pathways linked to OA pain (Supplementary Table 3C).However, an important mediator of pain in human and animal osteoarthritis, Nerve Growth Factor (NGF), although expressed in all our osteocyte cultures, was not regulated by load (Padj = 0.067).Expression of NGF has previously been reported to be upregulated in osteoblasts by physiological mechanical forces although it was not detected in osteocytes (94).In addition, several of the proposed mediators of OA pain such as neuropeptide Y (NPY), and substance P were not expressed in our model.However, Neuropeptide Y 1 receptor (NPY1R) was expressed and mechanically downregulated in our osteocyte model, Timeline of gene expression changes that occur following the differentiation of mature osteoblasts to mature osteocytes.Bars indicate cell stage (blue) showing the relative temporal expression of various osteogenic markers underneath (dark blue) during the transition from osteoblast to osteocyte with their respective protein role above (red).Adapted from ( 72) and (10).
consistent with the mouse osteocyte transcriptome (60) but we did not detect NPY itself, also in osteocyte transcriptome.NPY is released by osteocytes, is important in balancing adipocyte and osteoblast differentiation and mediates its effects via receptors, NPY1R and NPY2R which localise to pain centres in the nervous system (95).NPY influences nociceptive signalling in neuropathic and inflammatory pain (96).It is an important regulator of bone homeostasis and its expression is increased in osteoarthritic synovium and the concentration in synovial fluid from osteoarthritic patients positively correlates with pain scores, and is increased in late stages of OA (97).Since osteocyte NPY acts through NPY1R to suppress osteogenesis and promote adipogenesis of bone marrow stem cells, the effect of the mechanically induced 70% reduction in NPY1R expression that we observed requires further study.We and others have previously shown that glutamate receptors are involved in neural responses to inflammatory pain and that glutamate receptor antagonists alleviate pain and degeneration in animal models (13,22).AMPA (GRIA1/3/4), kainate (GRIK1/2/4/5) and NMDA (GRIN2A-D and 3A) ionotropic glutamate receptor subunits mRNAs were expressed in our osteocyte model.Of these GRIN2D and GRIK2 were identified in the osteocyte signature and involved with IDSWG3 maturation and GRIA3, GRIK5, GRIN2A and 3A in the osteocyte transcriptome (60).Kainate (GRIK2, 2-fold) and AMPA (GRIA3, 0.4-fold; GRIA4 20-fold) receptor subunits were mechanically regulated in osteocytes.This, along with reports of in vivo osteocyte expression of glutamate receptor proteins (AMPAR2 (13, 22); GRIN1, GRIA1/4 ( 14)), indicates a potential for mechanically-regulated glutamatergic signalling in osteocytes.Glutamate transporters (SLC1A1-3/EAAT1-3) were expressed by osteocytes consistent with previous reports in vivo (EAATs 1 and 2; SLC1A1 and 3 (60).Mechanical loading upregulated osteocyte SLC1A1/EAAT3 mRNA expression 3-fold but the mechanicallyinduced down regulation of SLC1A3/EAAT1 observed in osteocytes in vivo (16) was not recapitulated in this model.Upregulation of EAAT3 may explain why osteocytic glutamate release over 24 hours was reduced by mechanical loading.Increased concentrations of glutamate present in joint fluids in osteoarthritis, after joint injury, and at onset of joint inflammation (13, 22), are associated with pain and joint pathology (reviewed in (98)).The mechanically induced control of extracellular glutamate concentrations and regulation of glutamate receptor and transporter expression implicates osteocytes in the regulation of osteoarthritic pain and pathology although the effect of glutamate receptor activation in osteocytes is unknown.

Bone remodelling
Mechanical loading of osteocytes in our 3D model regulated 40% (23% upregulated; 16% downregulated) of the 217 genes known to be involved in gene ontology pathways linked to bone remodelling (Supplementary Table 3C).
The nearly 6-fold increase in OPG protein release after loading, without detectable changes in RANKL protein expression reveals a potential mechanically-induced inhibition of osteoclastogenesis and bone resorption consistent with the role of osteocytes in regulating bone remodelling in response to their mechanical environment (99).Mechanically regulated genes reflected pathways associated with RANKL signalling and abnormal bone cell function in a range of skeletal diseases, involving osteoclast activation and dysregulation of WNT signalling in osteoblasts.Interestingly similar pathways were identified as those that differed in the osteocyte transcriptome with age and sex (60).

Inflammation
Pathophysiological loading of osteocytes in our 3D model expressed 371 genes involved in acute and chronic inflammation with 27% of these genes being mechanically regulated (19% upregulated, 8% downregulated) (Supplementary Table 3D).NFkB signalling was particularly activated by mechanical load consistent with osteocytic loading playing a role in activating important proinflammatory and apoptotic pathways.In addition to the gene changes, several cytokine proteins were released by our osteocyte model, with mechanical loading reducing expression of both pro-resorptive (GM-CSF, IL-6, and RANTES) as well as antiinflammatory (MCP-1, IL-8, IP-10) proteins.This indicates that osteocytes can directly link mechanical loading to inflammation potentially mediating pathological processes in mechanically driven diseases such as osteoarthritis.

IL-6 regulate readouts of osteoarthritis in our 3D model of osteocytes
IL-6 and its soluble receptor were used at concentrations reported in human synovial fluids from patients with OA (32) or following knee injury (31) to stimulate inflammation in our 3D osteocyte model.Although IL6/sIL6r downregulated the AMPA glutamate receptor mRNA, GRIA1, treatment with IL-6 had limited effects on the glutamate signalling pathway and did not regulate expression of either glutamate transporters or glutamate release.Conversely, IL-6 treatment did affect markers of bone remodelling, halving OPG expression and reducing indicators of bone formation (BGLAP) and mineralisation (ALPL) indicative of increased resorption and reduced osteogenesis.Treatment with IL-6/sIL-6r resulted in the release of both pro-and anti-inflammatory cytokines causing a >10-fold increase in IL-12p70 and IL-4 and approximately doubling TNF-a and IFN-g protein release by osteocytes in our model all of which have been implicated in OA [reviewed in (100)] (30, 32),].IL-4 and IL12 have been shown to be antiosteoclastogenic and exert anti-resorptive effects on bone (101, 102).These data show that osteocytes in our 3D model respond to an inflammatory stimulus to modulate readouts of OA including markers of pain, bone remodelling and inflammation.

Links to human OA
We analysed our mechanosome dataset to determine whether genes that have been specifically identified as effector genes in OA patients were differentially expressed.A genome-wide association study (GWAS) meta-analysis across 13 international cohorts (826,690 individuals, 177,517 with osteoarthritis) identified 77 putative effector genes by analysis of functional genomics, mapping, eQTL, and associations with animal and human musculoskeletal and neuronal phenotypes (65).Of the 77 genes with strong evidence as effector genes (score3+), 83% were expressed in our osteocyte model with nearly half mechanically regulated (17 upregulated; 15 downregulated) (Supplementary Table 6A).The druggable genome database (reference ( 103) in ( 65)), revealed that twenty tier 1 (approved/clinical-phase drugs), five tier 2 (binding partners to approved drug targets) and twenty tier 3 (druggable pathways) were associated with these 77 putative causal genes.Of the 32 mechanically regulated genes in our model, identified as 'effector genes', 10 represented potential druggable targets (tier 1: TGFB1, TNC, CTSK, NOS3; tier 3A: GDF5, LTBP3, SERPINF1, NOG, LTBP1; tier 3B: MGP).A GWAS study by Tachmazidou et al. found 9 genes underlying monogenic forms of bone development diseases and ten likely early OA effector genes (104).Six genes underlying monogenic diseases were expressed in our dataset (2 upregulated by load) and 7/10 OA effector genes were present (2 upregulated and 2 down regulated by load) (Supplementary Table 6B).In addition, a search of the OMIM ® database (https://omim.org/) using the terms 'osteoarthritis AND bone revealed 97 gene entries with a phenotype related to OA; 26 of the associated genes were regulated by load in our dataset, 16 up regulated and 10 down regulated and another further 38 genes were present but not regulated (Supplementary Table 6C).A recent study by Zhou et al. (105) also reported links between OA and mechanical responsive osteocyte genes including POSTN, NID2, and ASPN; all regulated by load in our dataset.Collectively, this data shows that mechanical loading of osteocytes in our model regulates the expression of several genes shown to be important in human osteoarthritis susceptibility and potential treatment.

Limitations
There are several limitations to our model.The type I collagen gel is not mineralised or organised as it would be in vivo and represents newly formed osteoid; with time this may become change and future studies could examine this.This means it does not exhibit the physical and chemical properties of bone and whilst the osteocyte-like cells appear to have good molecular and morphological phenotypes, these will not be completely the same as those in vivo.We have applied compressive strain based on measurements of the gel under loading, and mimicked strains observed in rodent bones (45).The mechanical environment experienced by the cells is not the same as it would be in a mineralised bone since the osteocytes are not located within lacunae with fluid flowing through them under load and the strain on the cell and the way the strain changes over the culture period is not defined.The cells are subjected to compressive, tensile stretch and fluid low all of which will influence the cell's responses; further work, beyond the scope of the current study, could perform finite element modelling to help clarify this.Finally, the load was a single load of 3000 cycles that occurred over approximately 5 minutes.This could reflect joint trauma, or high strains caused by an episode of abnormal loading through the joint.The intention is to identify osteocyte derived mechanoresponsive signals that could contribute to disease processes in osteoarthritis, not to model the chronic disease and all joint tissues.We only used osteocytes in the current study; future work could co-culture osteoblasts and/or osteoclasts to mimic more closely the bone environment or other cell types to investigate tissue interactions.PCA analysis revealed a potential variation in differentiation status across cultures, with the gene expression profile of sample C2 varying somewhat from C3-5 and the L2 and L3 response to load varying from L4-6.This variation is likely due to subtle differences in cell density or distribution when embedding within the type I collagen, leading to a delay in cell-cell interactions, network connectivity and inhibition of cell division (10).Despite this, we have shown that our cells express an osteocyte phenotype using several approaches including RT-qPCR expression of osteocyte markers, high homology to the osteocyte signature, which are genes enriched in osteocytes relative to bone marrow and other osteoblast lineage cells (60), and protein expression of the mechanosensing osteocyte protein, sclerostin; the cells would not express high levels of sclerostin.A timeline of gene expression changes that occur following the differentiation of mature osteoblasts to mature osteocytes is shown in Figure 7 [adapted from (10,72)].

Conclusion
We have developed a reproducible model of human osteocyte like cells that express >90% of the genes in the osteocyte transcriptome signature.Mechanical loading and inflammatory stimulation regulated many genes and proteins implicated in osteoarthritis symptoms of pain as well as inflammation and degeneration underlying disease progression.Nearly half of the genes classified as 'effectors' in GWAS were mechanically regulated in this model.This model reveals that osteocyte mechanobiology plays an important role in osteoarthritic pathology.The model will be useful in identifying new mechanisms underlying bone and joint pathologies and testing drugs targeting those mechanisms.

1
FIGURE 1 Osteocyte-like cells differentiated from Y201 mesenchymal stem cells are viable with an appropriate phenotype.Y201 mesenchymal stem cells were expanded in culture (A) prior to embedding in type I collagen gels (B-D).Within 1 hour of culture (B), cells were beginning to send out dendritic processes which were evident by 4 hours (C).A 3D network of interconnecting cells was clearly evident by day 7 (D).Phalloidin staining of cells revealed dendritic processes (E, F); green) and sclerostin expression (G, H); green) at day 7 of culture.DAPI nucleus = blue.Images were captured using a x4 (A), x10 (E), x20 (B-D, G, H), and x40 (F) objective.Scale bar = 20µm B-D, G,H; 100 µm (A, E, F).

3
FIGURE 3 Regulation of markers of pain and bone remodelling markers in the 3D osteocyte model.The amount of (A) glutamate and (B) OPG released into the media from loaded osteocytes at 1 hour and 24 hours post-load (n=3/treatment, GLM ANOVA and Tukey's posthoc tests).