Hyper-physiologic mechanical cues, as an osteoarthritis disease relevant environmental perturbation, cause a critical shift in set-points of methylation at transcriptionally active CpG sites in neo-cartilage organoids

Background: Osteoarthritis (OA) is a complex, age-related multifactorial degenerative disease of diarthrodial joints marked by impaired mobility, joint stiffness, pain, and a significant decrease in quality of life. Among other risk factors, such as genetics and age, hyper-physiological mechanical cues are known to play a critical role in the onset and progression of the disease (1). It has been shown that post-mitotic cells, such as articular chondrocytes, heavily rely on methylation at CpG sites to adapt to environmental cues and maintain phenotypic plasticity. However, these long-lasting adaptations may eventually have a negative impact on cellular performance. We hypothesize that hyper-physiologic mechanical loading leads to the accumulation of altered epigenetic markers in articular chondrocytes, resulting in a loss of the tightly regulated balance of gene expression that leads to a dysregulated state characteristic of the OA disease state. Results: We showed that hyper-physiological loading evokes consistent changes in ML-tCpGs associated with expression changes in ITGA5, CAV1, and CD44, among other genes, which together act in pathways such as anatomical structure morphogenesis (GO:0009653) and response to wound healing (GO:0042060). Moreover, by comparing the ML-tCpGs and their associated pathways to tCpGs in OA pathophysiology, we observed a modest but particular interconnected overlap with notable genes such as CD44 and ITGA5. These genes could indeed represent lasting detrimental changes to the phenotypic state of chondrocytes due to mechanical perturbations that occurred earlier in life. The latter is further suggested by the association between methylation levels of ML-tCpGs mapped to CD44 and OA severity. Conclusion: Our findings confirm that hyper-physiological mechanical cues evoke changes to the methylome-wide landscape of chondrocytes, concomitant with detrimental changes in positional gene expression levels (ML-tCpGs). Since CAV1, ITGA5, and CD44 are subject to such changes and are central and overlapping with OA-tCPGs of primary chondrocytes, we propose that accumulation of hyper-physiological mechanical cues can evoke long-lasting, detrimental changes in set points of gene expression that influence the phenotypic healthy state of chondrocytes. Future studies are necessary to confirm this hypothesis.


Introduction
Osteoarthritis (OA) is a complex, age-related multifactorial degenerative disease of the diarthrodial joints marked by impaired mobility, joint stiffness, pain, and a signi cant decrease in quality of life.Among other risk factors, such as genetics and age, hyper-physiological mechanical cues are known to play a critical role in the onset and progression of the disease (1).OA is characterized by an imbalance in the articular chondrocytes' anabolic and catabolic activities, impacting the integrity of the cartilage.
Previous studies have characterized the deregulated signaling pathways in articular chondrocytes in response to hyper-physiological mechanical cues with transcriptome-wide differential expression analyses.These studies show that injurious mechanical cues signi cantly enhance cell apoptosis (2) and cellular senescence (3), increase catabolic gene expression (4), and reduce matrix production (5), whereas physiologic mechanical loading induces a broad anabolic response in the transcriptome that is associated with increased matrix formation (6).Post-mitotic cells, such as articular chondrocytes, heavily rely on methylation at CpG sites to adapt to environmental cues and maintain phenotypic plasticity (7).However, these long-lasting adaptations may eventually have a negative impact on cellular performance (8,9).We hypothesize that hyper-physiologic mechanical loading leads to the accumulation of altered epigenetic markers in articular chondrocytes resulting in a loss of the tightly regulated balance of gene expression to a dysregulated state characteristic of the OA disease state.
Here, we aimed to study the effect of hyper-physiological mechanical stress on changes in DNA methylation-driven set points of epigenetically regulated gene expression that potentially contribute to OA-related loss of the chondrocytes' epigenetic controlled healthy maturational arrested phenotypic state.
To this end, we employed an established human induced pluripotent stem cell (hiPSC)-derived cartilage organoid model and studied the methylome and transcriptome-wide changes in response to previously assessed hyper-physiological mechanical loading conditions (10).Using these techniques, we show that changes in the epigenetic set point of transcription in chondrocytes responding to hyper-physiological loading overlap with OA pathophysiology, further underlining their mutual role in evoking aberrant chondrocyte cellular functions.

Results
Characterization of experimental set-up to test the epigenome-and transcriptome-wide effects of hyperphysiological loading.
To test the effects of hyper-physiological loading on the epigenetically regulated transcriptome, we employed a human induced pluripotent stem cell (hiPSC)-derived neo-cartilage organoid model.hiPSCs were differentiated into chondrocytes using a previously established chondrogenic differentiation protocol (11).To study the response of chondrocytes to hyper-physiological mechanical loading conditions we have applied two different models that are commonly used in osteoarthritis research: (1) chondrocyte derived (spherical) neo-cartilage constructs, and (2) chondrocytes embedded in (cylindrical) agarose constructs.To distill the most consistent effects we employed both models and performed a meta-analysis on the methylome-wide changes in response to hyper-physiological mechanical loading.
Successful differentiation towards chondrocytes and the production of neo-cartilage was con rmed by protein immunolabeling of collagen II (COLII) and collagen VI (COLVI), as well as staining for sulfated glycosaminoglycans (sGAGs) (Fig. S1A).
To test the hypothesis that injurious mechanical stress can alter the epigenetically regulated transcriptome, both of the organoid models were exposed to hyper-physiologic mechanical loading (n = 26) alongside unloaded controls (n = 27) (10).Methylome-and transcriptome-wide pro les, together with real-time quantitative polymerase chain reaction (RT-qPCR), immunohistochemistry (IHC), and dimethyl methylene blue (DMMB) assays were measured 12 hours after mechanical stimulation.First, we characterized the response of neo-cartilage organoids to hyper-physiological loading conditions by targeted analysis using RT-qPCR of catabolic and anabolic cartilage markers and mechano-sensors (Table 1).Similar to other mechanically induced injurious in vitro and in vivo models of OA (3,12), expression of anabolic ADAMTS5 signi cantly increased in response to hyper-physiological loading.This nding suggests that hyper-physiological loading conditions induced a catabolic response in neocartilage organoids.Additionally, there was an increase in PIEZO1 expression, which is hypothesized one of the main transducers of hyper-physiological mechanical loading (13).Nonetheless, staining intensity of COLII and COLVI, as well as sGAG deposition normalized to DNA content, showed no signi cant change in response to hyper-physiological loading conditions (Fig. S1B-C).1A; Table S1) and plotted them across the genome annotated with the GENCODE basic V12 database (Fig. 1B) (14).Notable examples of highly signi cant DM CpG sites are cg27310485 annotated to calcium binding protein S100A2 (beta = .01,FDR = 7.8X10 − 4 ), cg16217885 annotated to in ammatory receptor IL1R1 (beta = 0.01, FDR = 1.2x10 − 3 ), cg12795959 annotated to SIRPB1 (beta = -0.03,FDR = 2.4x10 − 7 ).As shown in Fig. 1B, we recognized some skyscrapers suggesting differentially methylated regions (DMRs).Upon de ning DMRs as 3 or more DM CpGs with an inter-CpG distance of < 1kb and allowing for 3 non-DM CpGs in the complete DMR (15), these DMRs did not reach statistical signi cance.To gain insight into the functional aspects of DM CpGs, their enrichment to chromatin states was analyzed (16).As shown in Fig. 1C, DM CpG sites with hyper-physiological loading were signi cantly enriched within chromatin regions associated with active transcription start site (TSS) proximal promoter states TssAFlnk (OR = 0.87, FDR = 1.8X10 (beta=-0.0632,FDR = 8.8X10 − 5 ), cg13727613 annotated RUNX2 encoding a well-known transcription factor detrimental to cartilage, (beta=-0.021,FDR = 6.2X10 − 4 ), and cg23194024 annotated MSL1 encoding a component of the histone acetyltransferase complex responsible for the majority of histone H4 acetylation (beta=-0.032,FDR = 6.9X10 − 5 ) (Fig. 1D).

Transcriptionally active CpG-sites
To biologically interpret the DM CpG sites in response to hyper-physiological loading more speci cally, we next integrated a previously assessed RNA sequencing dataset (17) of the same experiment.To prioritize DM CpG that likely affect gene expression, we rst prioritized, among the DM CpG sites, those that mapped to genes within 200 or 1500 bp of the transcription start site (TSS200, TSS1500), located within the 3' or 5' UTR regions, or CpG sites that were exon bound.This resulted in 2492 CpG-gene pairs (Fig. 2; Table S2).Next, among the 2492 DM CpG sites, we prioritized those that mapped to a gene that was also differentially expressed in response to hyper-physiological loading conditions based on FDR correction for the number of genes overlapping with mapped CpGs.This selection of the most likely transcriptionally active DM CpGs, henceforth de ned as mechanical loading induced, transcriptionally active CpGs (ML-tCpGs), consisted of a total of 208 ML-tCpGs that mapped at TSS200 (15.4%),TSS1500 (28.5%), 3'UTR (13.65%), 5'UTR (32.2%), and 1st Exon (5.78%) or were Exonbound (4.3%) (Table S3).As shown in Fig. 2 and Table S3, these 208 ML-tCpGs were connected to 169 unique differentially expressed (DE) genes.As shown in Fig. 3A, 57 of the 169 genes showed a strong protein-protein interaction (FDR = 2.4x10 − 5 ) as determined by STRING-DB.Examples of highly connected ML-tCpGs-genes within this network are HSPA1B, CD44, and CAV1 (Fig. 3B-C) To gain insight into the biological processes that are affected by the ML-tCpGs responding to hyper-physiological loading conditions, we performed pathway enrichment analysis of gene ontology biological processes (GO BP), KEGG, and Reactome (Table S4) on all 169 unique genes that were associated with a ML-tCpG.Pathways such as anatomical structure morphogenesis (e.g., ANGPTL4, WNT9A, HMGA2), wound healing (e.g., HSPB1, SERPINE2, FERMT1), and caveola assembly (e.g., CAV1, CAV2, PACSIN2) were enriched.
Overlapping key nodes in the network of t-CpGs affected by mechanical loading and OA pathophysiology To explore the role of the identi ed ML-tCpG-gene pairs in OA pathophysiology, we examined the overlap of our ML-tCpG-gene pairs with those previously reported tCpG-gene pairs associated to OA pathophysiology, i.e. differentially expressed between lesioned and preserved cartilage from OA patients who underwent a joint replacement surgery (18).Although, the overlap was modest (n = 8 of 142 OA-tCpG-gene pairs) (Fig. S3), the ML-tCpG genes that did overlap with OA-tCpGs, CD44, ITGA5, and CAV1, are central and particularly interconnected genes in the network of both ML-tCpG-gene pairs and OA-tCpGgene pairs (Fig. 3A; Fig. S2).Taken together, we showed that hyper-physiological loading resulted in changes in tCpG-gene pairs that are both located within gene regulatory chromatin states and that are central and responsive to changes that occur during OA pathophysiology.

Correlation of tCpGs with clinical OA phenotypes
Finally, we set out to gain insight into the clinical relevance of altered epigenetic control of the identi ed ML-tCpG.Therefore, we determined the association of ML-tCpGs overlapping with OA-related epigenetically regulated genes (7) to phenotypic traits in the respective preoperative radiograph of the OA patient (the RAAK study).This subset of ML-tCpGs was regulating PFKP, CD44, CAV1, SVIL, and CY1BP1.Next, to assess how methylation levels of these genes relate to OA severity, methylation levels in preserved cartilage from joint replacement surgeries of these overlapping ML-tCpGs were correlated with OA severity as determined by the Kellgren-Lawrence grading scale (KL-score), adjusted for BMI and age.
We found that both CD44 (beta=-0.053,P = 0.048) and PFKP (beta=-0.111,P = 3.12X10 − 3 ) methylation levels were negatively associated with the KL-score, suggesting indeed that mechanical loading-induced alterations in epigenetic set points of expression are associated to the OA disease state of the patient.

Discussion
The goal of this study was to determine the effect of hyper-physiological mechanical loading on changes in stable set points of epigenetically regulated gene expression (ML-tCpGs) that could contribute to the long-lasting, detrimental changes in chondrocytes that are characteristic of an OA phenotype.To this end, we employed two human induced pluripotent stem cell (hiPSC)-derived neo-cartilage organoid models for robust readouts and studied the methylome-and transcriptome-wide changes in response to hyperphysiological mechanical loading conditions.We showed that hyper-physiological loading evokes consistent changes in ML-tCpGs associated with expression changes in ITGA5, CAV1, and CD44, among other genes, which together act in pathways such as anatomical structure morphogenesis (GO:0009653) and response to wound healing (GO:0042060).Moreover, by comparing the ML-tCpGs and their associated pathways to tCpGs in OA pathophysiology, we observed a modest but particular interconnected overlap with notable genes such as CD44 and ITGA5.These genes could indeed represent lasting detrimental changes to the phenotypic state of chondrocytes due to mechanical perturbations that occurred earlier in life.The latter is further suggested by the association between methylation levels of ML-tCpGs mapped to CD44 and OA severity.
Since injurious mechanical loading is considered an important driver of the onset and progression of the OA, here, we studied for the rst time whether injurious mechanical loading evokes stable, detrimental changes to chondrocyte phenotypic states.In doing so, we revealed that DM CpGs are particularly enriched in transcription start sites and enhancers suggesting that the DM CpG sites evoke changes in gene transcription.However, we cannot exclude the epigenetic regulatory effects of in trans t-CPGs.To allow biological interpretation of the DM CpGs in response to hyper-physiological loading, we integrated RNA sequencing data from the same experiment.We showed that signi cant differential CpG-gene pairs with hyper-physiological loading occurred, particularly near genes OA relevant genes (e.g., ITGA5, CD44, CAV1, WNT9A, and HMGA2).The relevance of mechanically induced expression of ITGA5 is that the binding of matrix fragments such as bronectin to integrin α5β1 heterodimer activates a pro-catabolic response (19).Also, CD44 plays a role in matrix catabolism by degrading hyaluronic acid in articular cartilage (20), while catabolic stress has been shown to upregulate CAV1, coding for caveolin-1, which has been linked to chondrocyte senescence (21).Genes like WNT9A as well as HMGA2 are reported OA risk genes (22).Nonetheless, the design of our study does not justify a direct causal relationship between DM at ML-CpGs and differential gene expression.To further con rm a direct causal relationship between CpG-speci c methylation levels and gene expression, a CRISPR-Cas9-DNMT/TET1 guided manipulation of methylation levels of the identi ed CpG sites mapped to CD44, PFKP, SVIL, and CY1BP1, is warranted.
Here we have combined genome-wide methylation and RNAseq analysis of hiPSC-derived chondrocytes either in deposited (spherical) neo-cartilage or embedded in (cylindrical) agarose, which currently are the most commonly used models in osteoarthritis research (6, 10, 13, 23).The advantage of the spherical neo-cartilage model is that it contains an extracellular matrix deposited by the chondrocytes.Herein the response of the chondrocyte to the mechanical perturbation likely re ects changes in the chondrocytematrix interaction.Additionally, it allows for evaluating responses of mechanical loading on sGAGs and other matrix constituents by histology.On the other hand, the strain distribution on the spherical pellets, hence chondrocytes, is less equal and could have reduced power or introduce bias.The advantage of the cylindrical-shaped model, for that matter, allows for an equal strain distribution, hence a more precise relation between the applied stress and the response of the chondrocytes to the deformation.By performing a meta-analysis on the methylome-wide landscape of these two models we aimed to distill the most consistent and robust effects.
Founded by the notion that early environmental challenges could evoke long-lasting changes to methylation at tCpG sites, resulting in detrimental changes to set points of gene expression, we sought evidence that ML-tCpGs are associated with previously assessed, OA-related differences in methylation in articular cartilage (OA-tCpGs) (7).We showed that genes such as CAV1, CD44, and ITGA5 appeared to be identically changed, highly interconnected, and central to both the OA-tCpGs and ML-tCpGs networks.As such, it is tempting to suggest that mechanical injurious loading can indeed contribute to stable and detrimental changes to the phenotypic state of chondrocytes eventually making the chondrocyte prone to an OA phenotype.Nonetheless, this hypothesis needs to be further investigated by con rming that CD44 upregulation leads to detrimental downstream effects of chondrocytes towards an OA chondrocyte phenotype.
After subjecting the neo-cartilage organoids to a single episode of mechanical injurious loading, the epigenetic and transcriptomic pro les were captured 12 hours later.Although this challenge resulted in many changes in methylation and associated expression, we did not observe changes in matrix content.
The fact that histological changes were not visible is likely due to the relatively early time point and/or the intrinsic insensitivity of identifying changes in protein levels.Moreover, although changes in methylation at CpG sites are generally considered stable and long-lasting ( 9), the fact that we only measured methylation at 12 hours post-stimulus did not allow us to con rm the duration of the change in methylation, and future studies may wish to address this question directly.

Conclusion
Together, the current study con rms that hyper-physiological mechanical cues evoke changes to the methylome-wide landscape of chondrocytes, concomitant with detrimental changes in positional gene expression levels (ML-tCpGs).Since CAV1, ITGA5, and CD44 are subject to such changes and are central and overlapping with OA-tCPGs of autologous cartilage, we advocate that accumulation of hyperphysiological mechanical cues can evoke long-lasting, detrimental changes in set points of gene expression that eventually affect the phenotypic healthy state of chondrocytes.Future studies are necessary to con rm this hypothesis.

Experimental design
The objective of the current study was to study the effects of hyper-physiologic mechanical loading conditions on the epigenetically regulated transcriptome.Here, we employed an hiPSC-derived neocartilage organoid model that was exposed to hyper-physiological mechanical loading conditions.These samples were then analyzed using 850k EPIC array and RNA sequencing.
hiPSC line and cell culture An hiPSC line as described earlier was used (11).In short, the RVR-hIPSC line was retrovirally reprogrammed from BJ broblasts and characterized.The hiPSCs were maintained under standard conditions (37 ℃, 5% CO 2 ) on Matrigel (Corning) coated plates and refreshed daily with TeSR-E8 medium (STEMCELL Technologies) upon reaching approximately 70% con uence.

Neo-cartilage organoid models
Two different chondrogenic constructs were used for downstream analysis; either these chondrogenic constructs were directly used for further experiments or they were dissociated using collagenase II, encapsulated in 2% w/v agarose at 30 million cells/ml, and cultured for 14 days with CD creating cylindrical shaped constructs.

Mechanical loading
The spherical shaped neo-cartilage constructs were mechanically loaded using a MACH-1 mechanical testing device (Biomomentum, Laval, Canada), at a rate of 5hz with 20% sinusoidal peak-to-peak strain for 10 minutes as described earlier (10).The cylindrical constructs were loaded with a custom-build mechanical loading device, with the same loading regime.

Histology and immunohistochemistry
Neo-cartilage samples were xed in 4% formaldehyde and embedded in para n.Sections were stained with Alcian Blue (Sigma-Aldrich) and Nuclear Fast Red (Sigma-Aldrich).Deposition of collagen VI and collagen II in the neo-cartilage constructs was visualized immunohistochemically using a polyclonal antibody for COL6A1 (abcam ab6588), a primary sub-unit of COLVI, and a polyclonal antibody for COL2A1 (abcam ab34712), a primary sub-unit of COLII., antigen retrieval was done by treating depara nized sections with proteinase K (5 µg/ml, Qiagen) and hyaluronidase (5 mg/ml, Sigma).Sections were incubated overnight with the primary antibodies, followed by incubation with a HRP conjugated secondary antibody (ImmunoLogic).Peroxidase binding for collagen VI was visualized using diaminobenzidine, and sections were counterstained with haematoxylin.

RT-qPCR
Per sample, two replicate neo-cartilage pellets were collected in TRIzol (Invitrogen™, Carlsbad, CA, USA), and RNA was isolated using the RNeasy Mini Kit (Qiagen, Venlo, the Netherlands) according to the manufacturer's protocol.DNA contamination was removed by treating the RNA with Rnase-Free DNase.RNA quality (A260/280: reactive (demonstrated to map to multiple places in the genome) were ltered out (30).The probes that were overlapping with rare SNPs (probes in TFBS that showed extreme methylation pattern) were ltered out (31).To minimize the unwanted variation within and between samples, we used the Functional Normalization method from the min .1.36.0 R package (32).We ran differential mean analysis using tmoderated statistics.Using the MEAL.1.20.3R package pipeline, which, relies on the lmFit from limma R package (design model=~ Loading).CpGs after Bonferroni correction P < 6.243109e-08 (0.05/800883) were considered signi cant.Strati ed analysis for each neo-cartilage construct was performed.These two datasets were then combined with a random effect meta-analysis using the metaVolcano R package The circus plot was produced using the Circlize 0.4.3R package (33).

Statistical analysis
For all data analysis except methylome data, we have used a generalized linear model including the factors hyper-physiological loading using R statistical software version 4.1.1.

Not applicable
Effect of hyper-physiological loading on the genome-wide methylation (A) A volcano plot of the methylome-wide response to hyper-physiological loading conditions.Blue dots denote CpG sites mapped to a gene body with increased methylation, FDR<0.01, and red dots represent CpG sites mapped to a gene body that are de-methylated in response to hyper-physiological loading conditions as determined by MEAL.(B) Manhattan plot of differentially methylated CpG sites with their genomic mapped genes.The horizontal red line represents the FDR<0.05threshold.(C) Enrichment of signi cant DMs within chromatin states; active transcription start site (TSS), proximal promoter states (TssA, TssAFlnk), a transcribed state at the 5′ and 3′ ends of genes showing both promoter and enhancer signatures (TxFlnk), actively transcribed states (Tx, TxWk), enhancer states (Enh, EnhG), and a state associated with zinc nger protein genes (ZNF/Rpts).The inactive states consist of constitutive heterochromatin (Het), bivalent regulatory states (TssBiv, BivFlnk, EnhBiv), repressed Polycomb states (ReprPC, ReprPCWk), and a quiescent state (Quies) (D) Notable examples of differentially methylated mapped CpGs in response to hyperphysiological mechanical loading conditions.The box plots represent the 25th, 50th, and 75th percentiles, and whiskers extend to 1.5 times the interquartile range.Individual samples are depicted by black dots in each graph.*FDR<0.05,**FDR<0.01,***FDR<0.001.

Table 1
Effects of mechanical loading on general OA markers (29)oAmp™ Optical 384-Well Reaction Plate (ThermoFisher Scienti c, Landsmeer, the Netherlands), using the QuantStudio™ Flex Real-Time PCR system (Applied Biosystems™, Foster City, CA, USA), with the following cycling conditions: 10 min 95 °C; 10 sec 95 °C, 30 sec 60 °C, 20 sec 72 °C (45 cycles); 1 min 65 °C and 15 sec 95 °C.Primer e ciency was tested using a cDNA dilution series, and primers were considered e cient with an e ciency between 90% and 110%.-ΔCtexpressionlevels were calculated using two housekeeping genes GAPDH and SDHA with the following formula: ΔCt=Ct (gene of interest) -Ct (average housekeeping genes).Both housekeeping genes were stably expressed in this model.Fold changes were calculated using the 2 -ΔΔCt method with ΔΔCt=ΔCt (MS) -ΔCt (Control).SOPs) based on the Illumina In nium II Protocol.To analyze methylation array data (MethylationEPIC 850k array), the MethylAid R script (28) with the default settings was used for data quality assessment.All samples showed a detected CpG above 95%.The min .v_1.36.0 R package(29)was used to pre-process the data.We removed any probe that have failed in one or more samples (p < 0.01).Probe level intensities were quantile normalized across samples prior to calculation of the ß-values.MethylToSNP was used to lter SNPs.This method looked for patterns in methylation array data and identi ed methylation probes with SNP-like patterns.The method removes outliers, which adds robustness to the analysis and is enabled by default.A con dence score was calculated to show how close the observed pattern of methylation beta values was to a canonical case of a SNP in a homozygously methylated CpG locus.Additionally, MethylToSNP can overlap the SNPs identi ed in methylation data with known SNPs from dbSNP.The probes that have shown to be cross-