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

Abstract

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.

1. 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 diagnosed with OA, the only available treatment options are joint replacement surgery and/or pain management. Therefore, there is a dire need for the discovery of biomarkers that can facilitate early detection of the disease and new therapeutic strategies for the prevention of PTOA. While many factors can influence the development of OA, injury-mediated OA holds the greatest promise for the development of effective pharmacologic interventions because a treatment can be administered at the time of surgery or immediately post-injury.

Acute joint trauma triggers several molecular events over the course of the first 1–2 weeks post-injury, which directly or indirectly contribute to the subsequent cartilage damage that is characteristic of OA. An understanding of these early molecular events provides a basis for identifying potential biologic targets for intervention to prevent subsequent joint degeneration [5]. Characterization of gene expression changes during OA development and progression at the whole genome level will provide novel mechanistic insights that could be translated into disease-modifying therapies. Numerous studies have used human biopsy samples to gain new insights about joint OA pathogenesis [6,7,8]; however, there are limitations in terms of the types of studies that can be conducted using human subjects. Human samples are usually obtained during knee replacement surgery and therefore represent late stages of the disease. To overcome this gap in knowledge, animal models allow us to investigate OA development longitudinally and are particularly well suited to studying early molecular events to derive new insights into the key factors contributing to disease progression. Using a noninvasive tibial compression (TC) injury model [9] that closely mimics anterior cruciate ligament (ACL) rupture in humans, we recently profiled the genes expression in knee joints from C57BL/6J mice during the onset and progression of PTOA and identified the molecular changes that characterize early and late stages of PTOA [10], including enhanced inflammatory responses at early timepoints and cartilage and bone remodeling at both early and late timepoints. We also noted that a majority of the transcriptional changes happened within the first few weeks post-injury.

In this study, we used the TC injury model to study early molecular events associated with PTOA development in three mouse strains with varying susceptibility to OA: STR/ort (highly susceptible), C57BL/6J (moderately susceptible), and MRL/MpJ (not susceptible). The STR/ort mice develop OA spontaneously early in life and show many human OA characteristics, including proteoglycan loss, extracellular matrix (ECM) degradation, and subchondral sclerosis [11]. They also exhibit osteophyte formation, a phenotype more readily seen in animal models where the joint is not stabilized. The MRL/MpJ is a mouse strain with exceptional abilities to heal wounds made in multiple tissues without the production of a fibrotic scar [12]. MRL/MpJs are protected from PTOA and do not develop degenerative joint changes following articular fracture [13]. It has been suggested that MRL/MpJ mice possess an intrinsic ability to regenerate articular cartilage, yet the molecular mechanisms responsible for this phenotype have yet to be revealed [14]. To characterize genes that contribute to increased OA susceptibility in STR/ort and resistance to PTOA in MRL/MpJ and to understand the molecular changes associated with early stages of PTOA development in these mice, we profiled gene expression in the knee joints of MRL/MpJ, C57BL/6J, and STR/ort by RNA sequencing (RNA-seq) at 0-day, 1-day, 1-week and 2-weeks post-injury. Understanding the molecular and genetic basis of enhanced OA susceptibility in STR/ort and resistance to PTOA in MRL/MpJ will improve our understanding of PTOA pathogenesis and may highlight new treatment options for PTOA or identify biomarkers that track disease progression.

This study identified 944, 2330, and 2702 genes differentially regulated in MRL/MpJ, C57BL/6J, and STR/ort, respectively, in response to ACL injury, including 553 genes that were shared by all strains. We identified increased, persistent inflammation, elevated catabolic activity, and elevated apoptosis as significant contributors to PTOA development. This study also identified several genes that may contribute to enhanced healing and tissue regeneration, including B4galnt2 and Tpsab1. Furthermore, this study identified several potential OA biomarkers, including Mamdc2 and Pxdn.

2. Results

2.1. 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).

Figure 1. Anterior cruciate ligament (ACL) injury leads to post-traumatic osteoarthritis (PTOA) in C57BL/6J and STR/ort but not in MRL/MpJ mice. Knee joints were injured at 10-weeks of age. (A) Timeline for histology, µCT, and RNA-seq sample collection [0-day (0D), 1-day post-injury (1D), 1-week post-injury (1W), 2-weeks post-injury (2W), 6-weeks post-injury (6W), and 12-weeks post-injury (12W)]. (B) µCT representation of injured and uninjured joints at 6-weeks post-injury. Darker regions in the injured scans depict osteophytes. (C) Osteophyte volume quantification (dark regions in the injured scans in (B) at 6-weeks post-injury. (D) Epiphyseal trabecular bone volume ratio of the distal femur was quantified and analyzed between injured and uninjured joints at 6-weeks post-injury. (E) Histological assessment of uninjured and injured joints at 12-weeks post-injury using Safranin-O and Fast Green staining (5× magnification). (F) OARSI scoring of histological sections of injured and uninjured joints at 12-weeks post-injury. Scale bar = 1 mm. * p < 0.05, ** p < 0.01, *** p < 0.001, ns not significant.

2.2. Characterizing Genes Related to Osteoarthritis Susceptibility

To identify transcripts that correlate with OA risk, we profiled the knee joints of 10-weeks old uninjured MRL/MpJ, C57BL/6J and STR/ort mice and examined all pair-wise comparisons (Figure 2A,B, Table S1). Four hundred ninety-seven genes were up-regulated in STR/ort compared to both C57BL/6J and MRL/MpJ (Figure 2A). This included 33 ‘inflammatory response’ -related genes (Il1b, Il6, Ccl2, Ccl7, Cxcl1, Cxcl2, etc.) (Figure 2C) and 78 genes associated with ‘apoptotic process’ including Egr1, Id3, Cebpb, Fos, and Jun (Figure S1A). Genes up-regulated in STR/ort compared to the other two strains also showed enrichment for ontology terms ‘response to wounding’ (40 genes), ‘response to oxidative stress’ (24 genes including Ptgs2, Rhob, Duox1, Areg and Plk3), ‘response to hypoxia’ (19 genes), and ‘ossification’ (25 genes)’. Enriched signaling pathways associated with these genes include ‘TNF signaling pathway’ (15 genes), ‘IL17 signaling pathway’ (14 genes) and ‘IL6-mediated signaling events’ (8 genes) (Figure 2D). Genes down-regulated in STR/ort compared to the other two strains showed enrichment for biological processes such as ‘leukocyte activation’ (29 genes), ‘hemopoiesis’ (27 genes) and ‘immune response’ (39 genes including Tlr1, Tlr5, Itgam, Itgad, Cd3d and Cd8a).

Figure 2. Genes associated with osteoarthritis susceptibility. (A) Overlap between genes upregulated in STR/ort (STR) compared to C57BL/6J (BL6) and MRL/MpJ (MRL) mice and in BL6 compared to MRL at 0-day. Genes upregulated in STR compared to both BL6 and MRL are shown in black ovals and genes downregulated in MRL (upregulated in both STR and BL6 compared to MRL) are shown in red ovals. (B) Overlap between genes downregulated STR compared to BL6 and MRL and in BL6 compared to MRL at 0-day. Genes downregulated in STR compared to both BL6 and MRL are shown in black ovals and genes upregulated in MRL (downregulated in both STR and BL6 compared to MRL) are shown in green ovals. (C) Inflammatory-response-related genes highly expressed in STR compared to BL6 and MRL at 0-day. (D) TNF, IL17, and IL6 signaling pathway-associated genes with high expression in STR compared to both BL6 and MRL. (E) Interferon signaling pathway-associated genes with low expression in MRL compared to both STR and BL6. (F) Matrix metalloproteinases (MMPs) with low expression in MRL compared to both STR and BL6.

MRL/MpJ had the lowest expression values for the majority of ‘inflammatory response’-related genes (Figure 2C). Genes with lowest expression in MRL/MpJ compared to other two strains (839 genes) also showed enrichment for ontology terms including ‘extracellular matrix organization’ (69 genes), ‘collagen catabolic process’ (18 genes), ‘response to wounding’ (71 genes), ‘immune response’ (98 genes), and ‘cytokine production’ (50 genes). Expression of interferon signaling pathway-associated genes, such as Ifitm1, Ifit3, Ifi27, Ifi35, Trim2, and Trim5, and matrix metalloproteinases (MMPs), including Mmp2, Mmp3, Mmp14, and Mmp19, were also significantly lower in MRL/MpJ compared to the other two strains (Figure 2E,F). Several genes involved in wound healing, including Postn, Fn1, Dcn, Ctgf, Col3a1, Col5a1, Plat, and Hbegf also exhibited lowest expression levels in MRL/MpJ (Figure S1B). Four hundred sixty-nine genes were upregulated in MRL/MpJ compared to both STR/ort and C57BL/6J (Figure 2B). This included 21 genes associated with ‘T cell activation’ (Cd3d, Cd4, Jak3 etc.). Other enriched ontology terms associated with genes upregulated in MRL/MpJ included ‘lymphocyte aggregation’ (21 genes), ‘defense response’ (49 genes), ‘heme metabolic process’ (7 genes), and ‘cellular ion homeostasis’ (23 genes). Enriched pathways associated with the upregulated genes included ‘heme biosynthesis’ (4 genes) and ‘complement cascade’ (6 genes).

2.3. Early Molecular Changes Associated with PTOA Development in STR/ort, C57BL/6J, and MRL/MpJ Mice

To identify the early molecular changes associated with PTOA development in STR/ort, C57BL/6J, and MRL/MpJ mice, we examined gene expression changes in injured joints compared to uninjured contralateral joints at 1-day, 1-week, and 2-weeks post-injury. At 1-day post-injury, we identified 441 (326 up; 115 down), 1590 (986 up; 604 down), and 735 (359 up; 376 down) differentially expressed genes (>1.5 fold) in the injured joints of MRL/MpJ, C57BL/6J, and STR/ort, respectively, including 201 genes commonly changed in all three strains (Figure 3A–D, 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”.

Figure 3. Early molecular changes associated with PTOA development. (A) Genes upregulated and (B) downregulated at 1-day (1D), 1-week (1W), and 2-weeks (2W) post-injury in STR/ort (STR), C57BL/6J (BL6), and MRL/MpJ (MRL). Overlap between genes (C) upregulated and (D) downregulated at 1-day post-injury. Overlap between genes (E) upregulated and (F) downregulated at 1-week post-injury. Overlap between genes (G) upregulated and (H) downregulated at 2-weeks post-injury. Genes upregulated in all three strains at (I) 1-day, (J) 1-week, and (K) 2-weeks post-injury. Majority of these injury-induced genes showed highest expression in injured STR.

Figure 4. Expression profiles of inflammatory cytokines. (A) Selected inflammatory cytokines with significantly higher expression in injured STR/ort (STR) compared to C57BL/6J (BL6) and MRL/MpJ (MRL) at 1-day post-injury. (B) Selected injury-induced genes with significantly higher expression in injured MRL compared to the other two strains at 1-day post-injury. I: Injured; U: Uninjured; 0D: 0-day; 1D: 1-day; 1W: 1-week; 2W: 2-weeks.

[ Table omitted. See PDF. ]

At 1-week post-injury, 654 (617 up; 37 down), 1203 (958 up; 245 down), and 2048 (1194 up; 854 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.

Figure 5. Matrix degrading enzymes showed highest expression in injured STR/ort. Expression of selected metalloproteinases in injured STR/ort (STR), C57BL/6J (BL6), and MRL/MpJ joints.

[ Table omitted. See PDF. ]

2.4. 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.

Figure 6. Candidates genes associated with enhanced healing in MRL/MpJ. (A) Genes differentially expressed in injured and uninjured MRL/MpJ (MRL) knee joints compared to both STR/ort (STR) and C57BL/6J (BL6) knee joints at all timepoints examined. (B) Elevated IHC expression of B4galnt2 in injured and uninjured joints of MRLs. Scalebar = 10 µm.

[ Table omitted. See PDF. ]

3. 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 6- and 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.

4. Materials and Methods

4.1. 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).

4.2. 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.

4.3. 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.

4.4. 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].

4.5. 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.

4.6. 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].

5. 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.

Supplementary Materials

Supplementary materials can be found at https://www.mdpi.com/1422-0067/19/9/2657/s1.

Author Contributions

Study design: G.G.L.; Data acquisition: J.C.C., A.S., M.E.M., D.K.M., B.A.C., S.H., and A.N.E. Data analysis and interpretation: A.S., B.A.C., and G.G.L., A.S., B.A.C. and G.G.L. wrote the manuscript.

Funding

A.S., J.C.C., B.A.C. and G.G.L. were supported by DOD grant OR130220 and LLNL grant 16-ERD-007. B.A.C. was supported by NIH/NIAMS grant AR062603.

Acknowledgments

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

ACL Anterior cruciate ligament

ECM Extracellular matrix

MMP Metalloproteinase

OA Osteoarthritis

PTOA Post-traumatic osteoarthritis

RNA-seq RNA sequencing

μCT Microcomputed tomography

References

1. Hunter, D.J.; Schofield, D.; Callander, E. The individual and socioeconomic impact of osteoarthritis. Nat. Rev. Rheumatol. 2014, 10, 437–441.

2. Lohmander, L.S.; Englund, P.M.; Dahl, L.L.; Roos, E.M. The long-term consequence of anterior cruciate ligament and meniscus injuries: Osteoarthritis. Am. J. Sports Med. 2007, 35, 1756–1769.

3. Thomas, A.C.; Hubbard-Turner, T.; Wikstrom, E.A.; Palmieri-Smith, R.M. Epidemiology of Posttraumatic Osteoarthritis. J. Athl. Train 2017, 52, 491–496.

4. Hannan, M.T.; Felson, D.T.; Pincus, T. Analysis of the discordance between radiographic changes and knee pain in osteoarthritis of the knee. J. Rheumatol. 2000, 27, 1513–1517.

5. Anderson, D.D.; Chubinskaya, S.; Guilak, F.; Martin, J.A.; Oegema, T.R.; Olson, S.A.; Buckwalter, J.A. Post-traumatic osteoarthritis: Improved understanding and opportunities for early intervention. J. Orthop. Res. 2011, 29, 802–809.

6. Ramos, Y.F.; Den, H.W.; Bovée, J.V.; Bomer, N.; Van, D.B.R.; Lakenberg, N.; Keurentjes, J.C.; Goeman, J.J.; Slagboom, P.E.; Nelissen, R.G.H.H.; et al. Genes involved in the osteoarthritis process identified through genome wide expression analysis in articular cartilage; the RAAK study. PLoS ONE 2014, 9, e103056.

7. Lambert, C.; Dubuc, J.E.; Montell, E.; Vergés, J.; Munaut, C.; Noël, A.; Henrotin, Y. Gene expression pattern of cells from inflamed and normal areas of osteoarthritis synovial membrane. Arthritis. Rheumatol. 2014, 66, 960–968.

8. Klinger, P.; Beyer, C.; Ekici, A.B.; Carl, H.D.; Schett, G.; Swoboda, B.; Hennig, F.F.; Gelse, K. The Transient Chondrocyte Phenotype in Human Osteophytic Cartilage: A. Role of Pigment. Epithelium-Derived Factor? Cartilage 2013, 4, 249–255.

9. Christiansen, B.A.; Anderson, M.J.; Lee, C.A.; Williams, J.C.; Yik, J.H.; Haudenschild, D.R. Musculoskeletal changes following non-invasive knee injury using a novel mouse model of post-traumatic osteoarthritis. Osteoarthr. Cartilage 2012, 20, 773–782.

10. Chang, J.C.; Sebastian, A.; Murugesh, D.K.; Hatsell, S.; Economides, A.N.; Christiansen, B.A.; Loots, G.G. Global molecular changes in a tibial compression induced ACL rupture model of post-traumatic osteoarthritis. J. Orthop. Res. 2017, 35, 474–485.

11. Staines, K.A.; Poulet, B.; Wentworth, D.N.; Pitsillides, A.A. The STR/ort mouse model of spontaneous osteoarthritis—An update. Osteoarthr. Cartilage 2017, 25, 802–808.

12. Heydemann, A. The super super-healing MRL mouse strain. Front. Biol. 2012, 7, 522–538.

13. Ward, B.D.; Furman, B.D.; Huebner, J.L.; Kraus, V.B.; Guilak, F.; Olson, S.A. Absence of posttraumatic arthritis following intraarticular fracture in the MRL/MpJ mouse. Arthrit. Rheum. 2008, 58, 744–753.

14. Fitzgerald, J.; Rich, C.; Burkhardt, D.; Allen, J.; Herzka, A.S.; Little, C.B. Evidence for articular cartilage regeneration in MRL/MpJ mice. Osteoarthr. Cartilage 2008, 16, 1319–1326.

15. Osterberg, A.; Thiem, D.; Herlyn, P.; Mittlmeier, T.; Frerich, B.; Müller-Hilke, B. Subchondral bone sclerosis and cancellous bone loss following OA induction depend on the underlying bone phenotype. Jt. Bone Spine 2017, 84, 71–77.

16. Cheng, C.H.; Leferovich, J.; Zhang, X.M.; Bedelbaeva, K.; Gourevitch, D.; Hatcher, C.J.; Basson, C.T.; Heber-Katz, E.; Marx, K.A. Keratin gene expression profiles after digit amputation in C57BL/6 vs. regenerative MRL mice imply an early regenerative keratinocyte activated-like state. Physiol. Genom. 2013, 45, 409–421.

17. Xu, R.; Singhal, N.; Serinagaoglu, Y.; Chandrasekharan, K.; Joshi, M.; Bauer, J.A.; Janssen, P.M.L.; Martin, P.T. Deletion of Galgt2 (B4Galnt2) reduces muscle growth in response to acute injury and increases muscle inflammation and pathology in dystrophin-deficient mice. Am. J. Pathol. 2015, 185, 2668–2684.

18. Pasold, J.; Engelmann, R.; Keller, J.; Joost, S.; Marshall, R.P.; Frerich, B.; Müller-Hilke, B. High bone mass in the STR/ort mouse results from increased bone formation and impaired bone resorption and is associated with extramedullary hematopoiesis. J. Bone Min. Metab. 2013, 31, 71–81.

19. Sokolove, J.; Lepus, C.M. Role of inflammation in the pathogenesis of osteoarthritis: Latest findings and interpretations. Ther. Adv. Musculoskelet. Dis. 2013, 5, 77–94.

20. Kyostio-Moore, S.; Nambiar, B.; Hutto, E.; Ewing, P.J.; Piraino, S.; Berthelette, P.; Sookdeo, C.; Matthews, G.; Armentano, D. STR/ort mice, a model for spontaneous osteoarthritis, exhibit elevated levels of both local and systemic inflammatory markers. Comp. Med. 2011, 61, 346–355.

21. Lewis, J.S., Jr.; Furman, B.D.; Zeitler, E.; Huebner, J.L.; Kraus, V.B.; Guilak, F.; Olson, S.A. Genetic and cellular evidence of decreased inflammation associated with reduced incidence of posttraumatic arthritis in MRL/MpJ mice. Arthritis. Rheum. 2013, 65, 660–670.

22. Heber-Katz, E.; Naviaux, R.K. The MRL Mouse: A Model of Regeneration and Cancer. In Murine Models, Energy Balance, and Cancer; Berger, N.A., Ed.; Springer International Publishing: Cham, Switzerland, 2015; Volume 10, pp. 47–64.

23. Gourevitch, D.; Clark, L.; Chen, P.; Seitz, A.; Samulewicz, S.J.; Heber-Katz, E. Matrix metalloproteinase activity correlates with blastema formation in the regenerating MRL mouse ear hole model. Dev. Dyn. 2003, 226, 377–387.

24. Tucker, B.; Klassen, H.; Yang, L.; Chen, D.F.; Young, M.J. Elevated MMP Expression in the MRL Mouse Retina Creates a Permissive Environment for Retinal Regeneration. Investig. Ophthalmol. Vis. Sci. 2008, 49, 1686–1695.

25. Burrage, P.S.; Mix, K.S.; Brinckerhoff, C.E. Matrix metalloproteinases: Role in arthritis. Front. Biosci. 2006, 11, 529–543.

26. Troeberg, L.; Nagase, H. Proteases involved in cartilage matrix degradation in osteoarthritis. Biochim. Biophys. Acta 2012, 1824, 133–145.

27. Flannelly, J.; Chambers, M.G.; Dudhia, J.; Hembry, R.M.; Murphy, J.; Mason, R.M.; Bayliss, M.T. Metalloproteinase and tissue inhibitor of metalloproteinase expression in the murine STR/ort model of osteoarthritis. Osteoarthr. Cartilage 2002, 10, 722–733.

28. Almodovar-Garcia, K.; Kwon, M.; Samaras, S.E.; Davidson, J.M. ANKRD1 acts as a transcriptional repressor of MMP13 via the AP-1 site. Mol. Cell Biol. 2014, 34, 1500–1511.

29. Samaras, S.E.; Almodóvar-García, K.; Wu, N.; Yu, F.; Davidson, J.M. Global deletion of Ankrd1 results in a wound-healing phenotype associated with dermal fibroblast dysfunction. Am. J. Pathol. 2015, 185, 96–109.

30. Liu, X.H.; Bauman, W.A.; Cardozo, C. ANKRD1 modulates inflammatory responses in C2C12 myoblasts through feedback inhibition of NF-kappaB signaling activity. Biochem. Biophys. Res. Commun. 2015, 464, 208–213.

31. Steinberg, J.; Ritchie, G.R.S.; Roumeliotis, T.I.; Jayasuriya, R.L.; Clark, M.J.; Brooks, R.A.; Binch, A.L.A.; Shah, K.M.; Coyle, R.; Pardo, M.; et al. Integrative epigenomics, transcriptomics and proteomics of patient chondrocytes reveal genes and pathways involved in osteoarthritis. Sci. Rep. 2017, 7, 8935.

32. Bouxsein, M.L.; Boyd, S.K.; Christiansen, B.A.; Guldberg, R.E.; Jepsen, K.J.; Müller, R. Guidelines for assessment of bone microstructure in rodents using micro-computed tomography. J. Bone Min. Res. 2010, 25, 1468–1486.

33. Trapnell, C.; Pachter, L.; Salzberg, S.L. TopHat: Discovering splice junctions with RNA-Seq. Bioinformatics 2009, 25, 1105–1111.

34. Kim, D.; Pertea, G.; Trapnell, C.; Pimentel, H.; Kelley, R.; Salzberg, S.L. TopHat2: Accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013, 14, 36.

35. Liao, Y.; Smyth, G.K.; Shi, W. FeatureCounts: An efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 2014, 30, 923–930.

36. Risso, D.; Ngai, J.; Speed, T.P.; Dudoit, S. Normalization of RNA-seq data using factor analysis of control genes or samples. Nat. Biotechnol. 2014, 32, 896–902.

37. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26, 139–140.

38. Chang, J.C.; Christiansen, B.A.; Murugesh, D.K.; Sebastian, A.; Hum, N.R.; Collette, N.M.; Hatsell, S.; Economides, A.N.; Blanchette, C.D.; Loots, G.G. SOST/Sclerostin Improves Post Traumatic Osteoarthritis and Inhibits MMP2/3 Expression After Injury. J. Bone Min. Res. 2018, 33, 1105–1113.

39. Chen, J.; Bardes, E.E.; Aronow, B.J.; Jegga, A.G. ToppGene Suite for gene list enrichment analysis and candidate gene prioritization. Nucleic Acids Res. 2009, 37, 305–311.

40. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 2003, 13, 2498–2504.