Comparative Transcriptomics Identifies Novel Genes and Pathways Involved in Post-Traumatic Osteoarthritis Development and Progression

Anterior cruciate ligament (ACL) injuries often result in post-traumatic osteoarthritis (PTOA). To better understand the molecular mechanisms behind PTOA development following ACL injury, we profiled ACL injury-induced transcriptional changes in knee joints of three mouse strains with varying susceptibility to OA: STR/ort (highly susceptible), C57BL/6J (moderately susceptible) and super-healer MRL/MpJ (not susceptible). Right knee joints of the mice were injured using a non-invasive tibial compression injury model and global gene expression was quantified before and at 1-day, 1-week, and 2-weeks post-injury using RNA-seq. Following injury, injured and uninjured joints of STR/ort and injured C57BL/6J joints displayed significant cartilage degeneration while MRL/MpJ had little cartilage damage. Gene expression analysis suggested that prolonged inflammation and elevated catabolic activity in STR/ort injured joints, compared to the other two strains may be responsible for the severe PTOA phenotype observed in this strain. MRL/MpJ had the lowest expression values for several inflammatory cytokines and catabolic enzymes activated in response to ACL injury. Furthermore, we identified several genes highly expressed in MRL/MpJ compared to the other two strains including B4galnt2 and Tpsab1 which may contribute to enhanced healing in the MRL/MpJ. Overall, this study has increased our knowledge of early molecular changes associated with PTOA development.


Introduction
Osteoarthritis (OA) is a painful degenerative joint disease that causes disability and diminishes the quality of life for millions of people worldwide [1]. Joint injury, particularly injuries to the anterior cruciate ligament (ACL), often result in post-traumatic osteoarthritis (PTOA) within 1-2 decades from the injury [2]. PTOA accounts for about 12% of all OA cases, yet the mechanisms contributing to PTOA after joint injury are not well understood, and there are currently no effective treatments available for PTOA [3]. Many people developing injury-or age-related OA do not show any symptoms until significant joint damage has occurred, and joint pain is not always indicative of OA [4]. For many

Evaluation of PTOA Development and Progression in STR/ort, C57BL/6J, and MRL/MpJ Mice Following ACL Injury
Injured and uninjured contralateral joints of STR/ort, C57BL/6J, and MRL/MpJ mice were phenotyped using histology and/or microcomputed tomography (µCT) at 6-and 12-weeks post-injury to assess tissue morphology. Gene expression was profiled by RNA-seq before injury and 1-day, 1-week, and 2-weeks post-injury ( Figure 1A). The joint damage was assessed by measuring the extent of osteophyte formation and the severity of cartilage degradation. Osteophyte formation was observed in all three strains by 6-weeks post-injury; MRL/MpJ injured joints had significantly less ectopic bone than the other strains ( Figure 1B,C). In all three strains, injured joints lost significant subchondral bone volume in the femoral epiphysis relative to the uninjured contralateral joints ( Figure 1D). Trabecular bone volume fraction (BV/TV) was significantly higher in STR/ort compared to the other two strains ( Figure 1D), and STR/ort had a significantly higher bone mass than the other two strains, consistent with prior publications [15]. At 12-weeks post-injury, STR/ort contralateral joints displayed significant proteoglycan loss and cleft down below the superficial and into the mid zone of the tibial cartilage, whereas MRL/MpJ and C57BL/6J contralateral joints had healthy cartilage ( Figure 1E). Injured joints of C57BL/6J and STR/ort exhibited severe cartilage erosion at 12-weeks post-injury ( Figure 1E,F). In contrast, the MRL/MpJ displayed an insignificant proteoglycan loss, suggesting that MRL/MpJ are protected from ACL injury induced cartilage damage. MRL/MpJ injured joints were significantly different than C57BL/6J and STR/ort injured joints, but no statistical difference was found between C57BL/6J and STR/ort injured joints ( Figure 1F).
to the other two strains ( Figure 1D), and STR/ort had a significantly higher bone mass than the other two strains, consistent with prior publications [15]. At 12-weeks post-injury, STR/ort contralateral joints displayed significant proteoglycan loss and cleft down below the superficial and into the mid zone of the tibial cartilage, whereas MRL/MpJ and C57BL/6J contralateral joints had healthy cartilage ( Figure 1E). Injured joints of C57BL/6J and STR/ort exhibited severe cartilage erosion at 12-weeks post-injury ( Figure 1E,F). In contrast, the MRL/MpJ displayed an insignificant proteoglycan loss, suggesting that MRL/MpJ are protected from ACL injury induced cartilage damage. MRL/MpJ injured joints were significantly different than C57BL/6J and STR/ort injured joints, but no statistical difference was found between C57BL/6J and STR/ort injured joints ( Figure 1F).    Table S2). Enriched biological processes associated with the upregulated genes included "extracellular matrix organization", "vasculature development", "cell migration", "angiogenesis", "response to wounding" and "inflammatory response" for all three strains and "granulocyte migration", "reactive oxygen species metabolic process", "cytokine secretion" and "response to tumor necrosis factor" for STR/ort alone. At 1-day post-injury, 59, 78, and 45 "inflammatory response"-related genes were upregulated in STR/ort, C57BL/6J, and MRL/MpJ, respectively. Several inflammatory cytokines, including Ccl2, Ccl7, Ccl8, Cxcl5, Il6, and Il33, were upregulated in all three genotypes in response to injury ( Table 1). The majority of these transcripts showed significantly higher expression levels in injured STR/ort compared to the other two strains ( Figure 4A). Furthermore, several immune/inflammatory response genes upregulated at 1-day remained elevated in STR/ort at later time points compared to uninjured controls; however, their expression returned to uninjured control level in C57BL/6J and MRL/MpJ strains by 1-2 weeks post-injury (Table 1). Only few injury-induced genes, including Ankrd1, Trpm1, and Fbxo32, showed highest expression in MRL/MpJ compared to the other two strains ( Figure 4B). Genes downregulated in STR/ort injured joints compared to uninjured controls showed enrichment for "biomineral tissue development", "muscle structure development" and "glycosaminoglycan metabolic process", whereas genes downregulated in MRL/MpJ injured joints compared to contralateral joints showed enrichment for "pyruvate metabolic process", "positive regulation of fatty acid oxidation", and "nucleotide phosphorylation".    down) genes were differentially regulated in injured joints of MRL/MpJ, C57BL/6J, and STR/ort, respectively, relative to respective uninjured controls ( Figure 3A,B,E,F). This included 269 genes commonly changed in all three genotypes (Table S2). At 2-weeks post-injury, 873 (778 up; 95 down) and 514 (501 up; 13 down) genes were differentially regulated in injured joints of STR/ort and C57BL/6J, respectively, whereas only 221 (215 up; 6 down) genes were found to be differentially regulated in injured joints of MRL/MpJ, and 183 of these genes overlapped with genes differentially expressed in the other two strains (Figure 3A,B,G,H). In addition, most of the genes differentially expressed at 2-weeks were upregulated, with less than 11% of differentially expressed genes being significantly downregulated in any strain (Table S2).
Although a large number of genes differentially expressed in response to injury were common to all three strains, STR/ort exhibited the highest expression values for a majority of these genes ( Figure 3I-K). Several genes associated with "extracellular matrix organization", "vasculature development", "response to wounding", "osteoblast differentiation", "ossification", and "collagen catabolic processes" were upregulated at 1-week and/or 2-weeks post-injury in all three strains. A number of catabolic enzymes, including Mmp2, Mmp3, Mmp19, Adamts1, and Adamts4, were upregulated in all three strains at 1-week and/or 2-weeks post-injury and the expression values of a majority of these catabolic enzymes were significantly higher in STR/ort compared to the other two strains (Table 2, Figure 5). At 1-week post-injury, genes associated with chondrocyte differentiation, including Sox9 and Runx1, were upregulated and several muscle-related genes, including Myh7, Myl2, Myl3, Myoc, Acta1, and Actc1, were downregulated exclusively in STR/ort and C57BL/6J (Table S2). Several members of Wnt signaling, a major signaling pathway involved in skeletal development and bone metabolism, including Wnt receptor Fzd2 and Wnt pathway inhibitors Sfrp1 and Sfrp2 were upregulated in all three strains at 1-week and 2-weeks post-injury (Table S2). It is likely that these genes play a significant role in cartilage and bone remodeling following ACL injury.

Potential Candidate Genes Associated with Enhanced Healing and Articular Cartilage Regeneration in MRL/MpJ
Compared to STR/ort and C57BL/6J, 204 genes were upregulated and 217 were downregulated in both injured and uninjured joints of MRL/MpJ at all timepoints examined (Table S3). ACL injury had little effect on the expression of majority of these genes ( Figure 6A). Using microarrays, Cheng et al. profiled genes involved in digit amputation response in MRL/MpJ and C57BL/6 at 0-day (pre-amputation), 3-days, 1-week, and 2-weeks post-amputation and identified genes differentially expressed in digits at various timepoints postamputation compared to 0-day in both strains as well as genes differentially expressed between strains [16]. To characterize the genes that may contribute to enhanced healing and/or regeneration in MRL/MpJ, we examined the overlap between genes differentially expressed in MRL/MpJ compared to the other strains in both knee joints and digits at 0-day, 1-week, and 2-weeks post-injury. Our analysis identified five genes upregulated in MRL/MpJ, including B4galnt2, Tpsab1, Vwa5a, and Aox4, and 16 genes downregulated in MRL/MpJ, including Mamdc2, Capg, Myoc, and Trim12a, compared to the other strains in both datasets at all timepoints examined ( Figure S2, Table 3). We further experimentally validated the differential expression of B4galnt2, a gene associated with muscle regeneration [17], in knee joints of all three strains. Immunohistochemistry confirmed that B4galnt2 ( Figure 6B) was highly expressed in both injured and uninjured joints of MRL/MpJ compared to other two strains.

Potential Candidate Genes Associated with Enhanced Healing and Articular Cartilage Regeneration in MRL/MpJ
Compared to STR/ort and C57BL/6J, 204 genes were upregulated and 217 were downregulated in both injured and uninjured joints of MRL/MpJ at all timepoints examined (Table S3). ACL injury had little effect on the expression of majority of these genes ( Figure 6A). Using microarrays, Cheng et al. profiled genes involved in digit amputation response in MRL/MpJ and C57BL/6 at 0-day (preamputation), 3-days, 1-week, and 2-weeks post-amputation and identified genes differentially     Coiled-coil domain containing 109B Thnsl2 Threonine synthase-like 2 Pccb propionyl Coenzyme A carboxylase, beta polypeptide Gpx3 glutathione peroxidase 3 Ezh1 enhancer of zeste 1 polycomb repressive complex 2 subunit Acsf2 Acyl-CoA synthetase family member 2 Pycard PYD and CARD domain containing

Discussion
STR/ort, C57BL/6J, and MRL/MpJ respond differently to knee joint injury. Here, we have introduced the first report that directly compares molecular and histological outcomes to a noninvasive ACL injury induced joint damage in these three strains. All three strains had deficits in epiphyseal trabecular bone in the injured joints and exhibited considerable osteophyte formation. STR/ort mice had degeneration in the contralateral joint and severe degeneration in the injured joint, whereas C57BL/6J mice had severe degeneration in the injured joint but not in the contralateral joint. By contrast, MRL/MpJ mice were almost completely protected from articular cartilage degeneration in this model. Consistent with previous reports [18], STR/ort showed higher trabecular bone volume fraction (BV/TV) compared to the other two strains. However, previous studies suggested that cartilage degeneration is independent of the underlying bone mass [15]-a hypothesis with diverging opinions in the field.
It has been suggested that the inflammatory response resulting from joint injury may be a significant factor in the progression of PTOA [19]. Studies have shown that STR/ort mice, a model for spontaneous osteoarthritis, exhibit elevated levels of both local and systemic inflammatory markers. Serum analysis showed elevated expression of several cytokines, including Il1β, Ccl4, and Il5, in STR/ort mice compared with that of CBA mice [20]. There is also evidence that MRL/MpJ mice have reduced inflammation, which may play a role in protecting these mice from PTOA [21]. Compared to C57BL/6 mice, MRL/MpJ mice had lower mRNA expression of Tnfα and Il1b in the synovial tissue and lower protein levels of Il1a and Il1b in the synovial fluid, serum, and joint tissues [21]. Consistent with these observations, our RNA-seq analysis showed that uninjured STR/ort joints expressed elevated inflammatory markers, including Il1b, Il6, Ccl2, and Cxcl1 and MRL/MpJ had the lowest expression values for these genes ( Figure 2C). Joint injury further amplified the expression of inflammatory-response-related genes at 1-day post-injury, which was greater in STR/ort than in the other two strains, and this inflammatory response persisted, while many of these genes reversed to pre-injury levels in C57BL/6J and MRL/MpJ shortly thereafter (Table 1, Figure 4A). This elevated, persistent inflammation may contribute significantly to the enhanced PTOA phenotype observed in STR/ort.
We observed that a number of genes associated with "T cell activation" were highly expressed in MRL/MpJ super-healer mice. However, immunodeficient MRL.RAG1 knockout mice were able to show complete ear hole closure, indicating that the regenerative response is not dependent on T or B cells in the ear [22]; it remains to be determined whether the same holds true for injured joints. We also observed higher expression of mast cell protease Tpsab1 in both injured and uninjured joints of MRL/MpJ compared to the other two strains at all timepoints ( Figure S2), which correlates elevated mast cells with the enhanced healing observed in this mouse strain. Another gene highly expressed in MRL/MpJ compared to the other two strains-B4galnt2-has been shown to play a role in skeletal muscle growth in response to acute muscle injury [17] (Figure 6B). We hypothesize that B4galnt2 also contributes to enhanced healing in MRL/MpJ, and future studies will address the role of this gene in PTOA. The regenerative healing characteristics of the MRL/MpJ strain can also be attributed to reduced expression of apoptosis-associated genes. STR/ort had the highest expression of several apoptosis-associated genes, including Fos, Jun, and Id3, whereas MRL/MpJ had the lowest expression ( Figure S1A).
In the mammalian adult, default response to injury involves inflammation, replacement of mature cells, and the formation of scar tissue. Healing in the MRL/MpJ appears more fetal-like with the formation of a blastema, healing without scarring, and with the replacement of lost tissue by functionally and architecturally normal tissue [23]. Remodeling and degradation of the ECM by MMPs is a key step in wound healing. Gourevitch et al. have shown that Mmp2 and Mmp9 protein levels were upregulated in the MRL/MpJ healing ear hole tissue compared with the C57BL/6 and that the MMP activity correlated with blastema formation in the regenerating ear holes [23]. It has also been shown that MRL/MpJ mice exhibit elevated levels of Mmp2, -9, and -14 in the retina compared to C57BL/6 and this elevated MMP expression creates a permissive environment for retinal regeneration in MRL/MpJ mouse [24]. Contrary to these prior findings, we determined that mRNA levels of Mmp2, Mmp3, and Mmp14 were significantly lower in the knee joints of MRL/MpJ compared to C57BL/6J and STR/ort joints ( Figure 5, Figure S3). Although the expression of these catabolic enzymes elevated in response to injury in all three strains, the expression of these genes remained significantly lower in the injured MRL/MpJ joints compared to injured C57BL/6J and STR/ort joints. These results may point out fundamental differences due to anatomical location and function. MMPs are critical for cartilage remodeling in joints, and elevated levels of these molecules in the joint are usually correlative of enhanced cartilage catabolic activity or degradation in OA [25,26]. Flannelly et al. have shown that mRNA levels of Mmp2, Mmp3, Mmp7, Mmp9, Mmp13, and Mmp14 were higher in STR/ort than in age-matched CBA mice at various ages [27]. Consistent with their findings, we found higher levels of MMPs and other matrix degrading enzymes, such as Adamts1, Adamts3, and Adamts4, in injured STR/ort relative to the other two injured strains ( Figure 5), suggesting that elevated levels of these matrix degrading enzymes may contribute to enhanced cartilage degradation in STR/ort, and lower levels in MRL/MpJ may be one mechanism by which this strain is resistant to PTOA.
Ankrd1, a transcriptional repressor of MMP13 [28], was highly upregulated in MRL/MpJ at 1-day post-injury, whereas only moderately changed in C57BL/6J and STR/ort ( Figure 4B). Global deletion of Ankrd1 resulted in delayed excisional wound closure [29]. Deletion of Ankrd1 also resulted in moderate downregulation of Mmp2 and Mmp14 [28]. Ankrd1 also plays an anti-inflammatory role through feedback inhibition of NF-κB transcriptional activity [30]. These findings suggest that Ankrd1 may play a role in protecting MRL/MpJ against injury-induced cartilage damage, possibly by keeping the MMP expression at a low level. MAM Domain Containing 2 (Mamdc2), a gene encoding an ECM protein, had extremely low expression in MRL/MpJ but was moderately expressed in the other two strains ( Figure S2). Mamdc2 was upregulated in injured STR/ort joints compared to uninjured contralateral joints at 2-weeks post-injury. Interestingly, we also found Mamdc2 to be upregulated in C57BL/6J at 6and 12-weeks post-injury [10]. Furthermore, MAMDC2 was significantly upregulated in human OA samples [31], which positions MAMDC2 as an ideal candidate biomarker for PTOA. Peroxidasin (Pxdn) is another potential OA biomarker identified in this study that was also upregulated in human OA [31]. Very little is known about the functions of these genes, and they warrant further investigation.
One limitation of our study is that we sequenced the whole joints instead of individual tissues of the joint, which makes it difficult to tease out the cellular source of the gene expressions observed. These challenges may be overcome by examining candidate proteins for their tissue-specific expression using other techniques such as immunohistochemistry. Another limitation is that we used the contralateral joint as controls instead of age-matched sham injured joints; this may have caused us to underestimate changes mediated by the injury that had systemic effects on both joints. Regardless, this study identified hundreds of genes and several new pathways that may contribute to PTOA pathogenesis and should be further evaluated in forthcoming studies. Our study provides novel insights into genes and molecular pathways involved in the early stages of PTOA development and identifies several putative candidate genes that may contribute to enhanced healing observed in MRL/MpJ. In addition, the data generated in this study could help facilitate future research in the identification and development of novel approaches to treat PTOA.

Animals and Tibial Compression (TC) Joint Injury
Right knee joints of 10-weeks-old STR/ort, C57BL/6J (Jackson Laboratory, Bar Harbor, ME, USA; Stock No: 000664), and MRL/MpJ (Jackson Laboratory, Bar Harbor, ME, USA; Stock No: 000486) mice were injured using a compressive load of 10-12N, as previously described [9,10]. Mice were anesthetized via isoflurane inhalation and placed in a prone position with right tibias vertically aligned between two platens for tibial compression. ACL rupture was produced via a single dynamic axial compressive load at 1 mm/s using an electromagnetic material testing machine (ElectroForce 3200, TA Instruments, New Castle, DE, USA). Buprenorphine analgesia was administered immediately post-injury (0.01 mg/kg). All animal experimental procedures were completed in accordance with the Institutional Animal Care and Use Committee (IACUC) guidance at Lawrence Livermore National Laboratory under approved protocol 168 (last approval date 5/21/2018).

Histological Assessment of Articular Cartilage and Joint Degeneration
Injured and uninjured (contralateral) joints were collected at 6-and 12-weeks (n > 5 per group) post-injury. Joints were dissected, fixed in 4% paraformaldehyde, decalcified using 0.5 M EDTA, infiltrated in increasing concentrations of isopropanol, equilibrated into mineral oil, and embedded into paraffin wax. Six-micrometer paraffin sections were stained on glass slides using 0.1% Safranin-O and 0.05% Fast Green using standard procedures (IHC World, Woodstock, MD, USA) and imaged using a Leica DM5000 microscope. Three blind reviewers independently assessed OA severity using a modified OARSI [scale to examine the medial compartment of injured and uninjured joints (sagittal views); grade scale 0-0.5 normal; 1-2 mild; 3-4 moderate; 5-6 severe cartilage damage]. All sagittal sections were collected on the medial compartment of the knee. For uninjured (contralateral) samples, the evaluation was performed by examining 6 µm serial sections between 400-600 µm in depth from the joint surface. For injured joints, due to the osteophyte and dislocation, histological evaluations were carried out on serial sections between 700-900 µm in depth. However, because MRL injured joints had reduced ectopic bone and joint dislocation, these samples were evaluated between 400-600 µm in depth to obtain comparable regions in the joint.

Microcomputed Tomography Analysis and Osteophyte Quantification
Whole knee joints (n > 5 per group) were scanned using µCT (SCANCO µCT 35, Brüttisellen, Switzerland) according to rodent bone structure analysis guidelines (X-ray tube potential = 55 kVp, intensity = 114 µA, 10 µm isotropic nominal voxel size, integration time = 900 ms) [32]. Trabecular bone in the distal femoral epiphysis was analyzed by manually drawing contours on 2D transverse slices. The distal femoral epiphysis was designated as the region of trabecular bone enclosed by the growth plate and subchondral cortical bone plate. Osteophyte volume in joints was quantified by drawing contours around all heterotopic mineralized tissue attached to the distal femur and proximal tibia as well as the entire patella, fabellae, and menisci; the patella, fabellae, and menisci of contralateral limbs were also contoured. Total mineralized osteophyte volume was then determined as the volumetric difference in mineralized tissue between injured and uninjured joints. Statistical analysis was performed using a paired t-test to compare injured and contralateral knees.

RNA Sequencing and Data Analysis
Injured and contralateral joints (n > 4 per group) were dissected and cut at the base edges of femoral and tibial joint regions with small traces of soft tissues to preserve the intact knee joint. The RNA was isolated and sequenced as previously described [10]. RNA-seq data quality was checked using FastQC (version 0.11.5) software. Sequence reads were aligned to the mouse reference genome (mm10) using TopHat (version 2.0.11) [33,34]. After read mapping, "featureCounts" from Rsubread package (version 1.22.2) [35] was used to perform summarization of reads mapped to RefSeq genes, and gene-wise read counts were generated. Genes were filtered from downstream analysis if they did not have read count of at least 2 in at least five libraries. RUVseq [36] was used to normalize data using 25 housekeeping genes (Supplementary Table S4). Differentially expressed genes were identified using edgeR (version 3.14.0) [37]. A gene was considered significantly differentially expressed when its false discovery rate (FDR) corrected p-value was less than 0.05 and fold change was greater than 1.5. Heatmaps were generated using heatmap.2 function in R package 'gplots'. Human OA RNA-seq data was obtained from Steinberg et al. [31].

Immunohistochemistry
Sagittal serial sections were stained utilizing primary antibodies against B4galnt2 (Novus Biologicals, Littleton, CO, USA). Trypsin/EDTA was used for antigen retrieval for 25 min at 37 • C. Antibody staining was performed as previously described [38]. Negative control slides were incubated with secondary antibody only. Stained slides were mounted with Prolong Gold with DAPI (Molecular Probes, Eugene, OR, USA). ImagePro Plus V7.0 Software and a QIClick CCD camera (QImaging, Surrey, BC, Canada) were used for imaging and photo editing.

Functional Annotation
Gene ontology analysis was performed using ToppGene [39] and enriched gene ontology terms and pathways (p-value < 0.01) were identified. Cytoscape was used for pathway visualization [40].

Conclusions
Our data suggest that prolonged inflammation and enhanced expression of matrix degrading enzymes may contribute to a severe PTOA phenotype. This study identified many new potential therapeutic targets, including B4galnt2, and potential OA biomarkers, including Mamdc2 and Pxdn. This study also highlighted several candidate genes that may contribute to enhanced healing and/or tissue regeneration.