Comparison of inbred mouse strains shows diverse phenotypic outcomes of intervertebral disc aging.

Abstract Intervertebral disc degeneration presents a wide spectrum of clinically degenerative disc phenotypes; however, the contribution of genetic background to the degenerative outcomes has not been established. We characterized the spinal phenotype of 3 mouse strains with varying cartilage‐regenerative potential at 6 and 23 months: C57BL/6, LG/J and SM/J. All strains showed different aging phenotypes. Importantly, LG/J mice showed an increased prevalence of dystrophic disc calcification in caudal discs with aging. Quantitative‐histological analyses of LG/J and SM/J caudal discs evidenced accelerated degeneration compared to BL6, with cellular disorganization and cell loss together with fibrosis of the NP, respectively. Along with the higher grades of disc degeneration, SM/J, at 6M, also differed the most in terms of NP gene expression compared to other strains. Moreover, although we found common DEGs between BL6 and LG/J aging, most of them were divergent between the strains. Noteworthy, the common DEGs altered in both LG/J and BL6 aging were associated with inflammatory processes, response to stress, cell differentiation, cell metabolism and cell division. Results suggested that disc calcification in LG/J resulted from a dystrophic calcification process likely aggravated by cell death, matrix remodelling, changes in calcium/phosphate homeostasis and cell transformation. Lastly, we report 7 distinct phenotypes of human disc degeneration based on transcriptomic profiles, that presented similar pathways and DEGs found in aging mouse strains. Together, our results suggest that disc aging and degeneration depends on the genetic background and involves changes in various molecular pathways, which might help to explain the diverse phenotypes seen during disc disease.


| INTRODUC TI ON
Despite the multifactorial aetiology of disc degeneration, one of the critical risk factors is aging (Chanchairujira et al., 2007;Cheung et al., 2009). This is becoming more important as increasing lifespan is contributing to a rise in age-dependent diseases. It is known that over the course of lifespan, discs experience decreased levels and quality of extracellular matrix proteins, compromised biomechanical properties, increase in inflammatory cytokine expression, catabolic processes and cell death (Caldeira et al., 2017;Ohnishi, Sudo, Tsujimoto, & Iwasaki, 2018). Similarly, senescence, herniations and calcifications have been correlated with degenerative phenotypes in the disc (Novais, Diekman, Shapiro, & Risbud, 2019). Accordingly, the term degenerative disc disease is an umbrella term that encompasses multiple phenotypes. It is therefore important to delineate the molecular basis of these distinct pathologies and their contribution to disc aging and degeneration.
The absence of an appropriate animal model is a barrier for a better understanding of the nuances of disc aging and degeneration.
Many experimental models, ranging from large to small animals, and from mechanically induced injury to genetic manipulations have been developed for this purpose (Caldeira et al., 2017;Choi et al., 2018;Novais et al., 2019;Yang, Leung, Luk, Chan, & Cheung, 2009). Injury models rely on acute trauma, usually by annular puncture with consequent inflammatory and immune responses. While injury models are of interest to study acute herniations and the role of inflammation in pathology, they do not recapitulate the chronic cases of human disc degeneration, especially those related to aging (Yang et al., 2009).
Likewise, genetic modification allows insights into how a single gene or protein contributes to the inflammation-driven disc herniations and/or degeneration but does not capture the multivariable complexity of the aging phenotype (Gorth, Ottone, Shapiro, & Risbud, 2019).
Nevertheless, models such as mice and rats are attractive for studying disc aging due to their genetic similarity with humans and the availability of tools for population genetics. Importantly, several studies have shown comparable changes in histology, mechanical properties, cellular responses and matrix remodelling during disc degeneration in mice and humans (Choi et al., 2018;Ohnishi et al., 2018). SM/J, a poor healer mouse strain, and LG/J, a super healer strain, have recently become an attractive system to study disc degeneration due to their opposing phenotypes (Choi et al., 2018) While SM/J mice showed early-onset, spontaneous disc degeneration, it is not known how their genetic background influences disc aging. To address these questions, we characterized the spinal phenotype of 6-month-old (6M) and 23-month-old (23M) SM/J, LG/J and C57BL/6 (BL6) mice, the most widely used strain in aging and disc research which exhibits poor healing capabilities compared to LG/J mice (Choi et al., 2018). Our results show that LG/J mice are prone to developing dystrophic disc calcification in caudal spine with aging and the aging phenotype is distinct between 3 strains. Transcriptomic analysis showed SM/J differs the most from the other strains at 6M, exhibiting highest grades of degeneration.
Interestingly, during aging, strains shared only a small subset of genes related to stress response, cell differentiation, inflammation-related pathways and cell metabolism with LG/J mice showing differences in pathways concerned with cell death, calcium homeostasis and endochondral processes. Signatures of calcium homeostasis, cell death and endochondral genes were also captured in hierarchical clustering analysis of transcriptomic data from human NP tissues. Together, the findings highlight the contribution of the genetic background to spectrum of degenerative phenotypes seen during aging and disc degeneration.

| LG/J mice have a higher prevalence of agedependent disc calcification than SM/J and BL6
To investigate how genetic background affects disc phenotype with aging, we analysed the spines of 6M and 23M BL6, LG/J and SM/J our results suggest that disc aging and degeneration depends on the genetic background and involves changes in various molecular pathways, which might help to explain the diverse phenotypes seen during disc disease.

K E Y W O R D S
Aging, calcification, intervertebral disc, LG/J, SM/J, transcriptome F I G U R E 1 LG/J mice present high prevalence of age-dependent disc calcification. (a-c) μCT showed higher prevalence and widespread distribution of disc calcification in caudal spine of 23M LG/J mice. Yellow, orange and red arrows (a) and colour scale (c) represent small, medium and large size calcifications, respectively. p ≤ .001***, χ 2 test, BL6 (n = 10), LG/J (n = 14) and SM/J (n = 8). (d-d') Disc calcifications in LG/J mice were present in both NP and AF. (e-e') Alizarin Red staining showed staining of the disc, not limited to the calcified nodules. (f) FTIR analysis showed comparable spectral profile between disc mineralization and vertebrae. Green arrow-amide peak -1,665/cm; red arrow-phosphate peak -960/cm; purple arrow-carbonate peak -870/cm; *phosphate peak-1,030/cm; ** tape artefact. (g-h') Phosphate and carbonate levels in the disc were comparable to the vertebrae. (i-i') Phosphate to protein ratio was higher in the disc. (j) Plasma measurements of calcium, phosphorous, TNAP, fetuin, glucose, creatine, sodium and albumin from 6M and 23M mice. Mann-Whitney test was used for comparing differences between the groups. ns = not significant; p ≤ .05*; p ≤ .01**; p ≤ .01**; BL6 (n = 10), LG/J (n = 14) and SM/J (n = 8). Scale bar d-e' = 200 µm mice by μCT and histology. μCT analysis showed high prevalence of disc calcification in caudal spine of 23M LG/J mice with a comparable mineral density to vertebral cortical bone (Figure 1a-c).
Calcification presented a wide distribution, size and morphology LG/J 6M LG/J 23M

Amide
Ca6-7, Ca7-8 and Ca8-9, whether this is due to differential loading at these levels remains to be seen. There was a lack of disc calcification in 23M-old BL6 and SM/J mice (Figure 1a, Figure S2a) and in LG/J lumbar spine ( Figure S2a). Alizarin Red staining showed presence of pronounced free calcium in 23M LG/J disc that was not restricted to the mineral nodules ( Figure 1d-e' and Figure S2b). FTIR-studies showed that this calcification was apatitic with a comparable spectral profile to the mineral in the vertebrae (Figure 1f). We analysed the peaks corresponding to phosphate-960/cm (Figure 1g  LG/J and SM/J mice showed increased activity of TNAP, without differences between the strains. Levels of fetuin-A were higher in LG/J but virtually undetectable in the SM/J mice. Plasma glucose was significantly decreased in SM/J mice as compared to BL6 and LG/J mice without any strain-dependent change in creatine, sodium and albumin (Figure 1j). These results suggested that there was no clear correlation between the plasma chemistry and strain-specific mineralization phenotype. We then decided to focus on analysing caudal disc aging due to this calcification phenotype and to avoid confounding contribution of axial loading on lumbar spine associated with variations in body sizes between three strains (Espinoza Orías, Malhotra, & Elliott, 2009).

| LG/J and SM/J mice have distinct phenotypes of age-related disc degeneration compared to BL6
We performed histological analysis of caudal disc of 6M and 23M BL6, LG/J and SM/J mice. There were differences in tissue archi-

| Aged inbred mice strains evidenced changes in extracellular matrix composition
We assessed the differences in disc matrix composition in 3 strains at . However, we also observed more pronounced peri-cellular staining of collagen X in the SM/J NP, in addition to higher AF levels compared to LG/J and BL6. These results showed that unique collagen signatures underscored disc aging phenotypes in these three strains (Choi et al., 2018).
We examined aggrecan content and functionality parameters.
LG/J and BL6 mice showed higher aggrecan abundance compared to LG

| Comparison of aging BL6 and LG/J mice uncovers biological pathways that are independent and dependent on the strain background
To explore the mechanisms characterizing disc aging and dystrophic mineralization in LG/J mice, we performed microarray analysis of NP tissues from 23M BL6 and LG/J mice and compared it to 6M data from respective strains. Lack of enough viable cells in SM/J NP precluded their analysis at 23M ( Figure S3) (Choi et al., 2018). The These results suggest that although the strains showed some commonality during aging, strain-dependent biological processes have a higher contribution to the transcriptomic signature of aged NP cells.

| LG/J shows increased cell death without changes in NP marker expression
We performed TUNEL assay to investigate whether there was in-
Levels of NP markers CA3 and GLUT1 were similar between the two strains but a complete absence of staining was seen in SM/J (Choi, 2018) ( Figure S7a).

| Inbred strain comparison shows differential effects of aging on vertebral bone parameters
Previous report suggested that LG/J mice have superior bone quality parameters compared to SM/J and BL6 mice (Sabsovich, 2008).
To assess the impact of aging on vertebral health, we compared the bone morphology in the caudal spine of 3 strains using μCT ( Figure   S1a). Although significant changes in the vertebral height were found between the strains, the disc height index (DHI) was comparable at 23M ( Figure S1b) showed lower values of trabecular bone than BL6 but thicker and more resistant cortical bone than the two strains.

| Degenerated human discs show distinct transcriptomic profiles and enriched biological pathways found in aging mouse strains
To understand the relationship between mouse models and human pathology, we re-analysed microarray data set of healthy and degenerated human NP tissues from a previous study by Kazezian et al. (2015) (GSE70362). Interestingly, the transcriptomic profiles of healthy versus degenerated discs did not cluster distinctly according to standard histological grades of degeneration ( Figure 6a). Therefore, using a hierarchical clustering of transcriptomes of all the available healthy (grade 1 and 2-in blue) and degenerated (grades 3, 4 and 5-in red) NP samples, with a cut-off of dissimilarity < 0.5 (Euclidean distance), we categorized the samples into 7 degenerative groups-red boxes, and 3 healthy groups-green boxes (Figure 6b). From these groups, we compared the cluster A (healthy cluster with most samples-5) to the remaining 7 degenerative clusters. This new analysis showed that each cluster organized together along the principal components (Figure 6c).
To further explore the differences and similarities of gene expression from these cluster comparisons, we performed enrichment biological processes analysis of both upregulated and downregulated DEGs. Due to the divergent transcriptome signature of each degenerative cluster, we found that while some enriched pathways were similar, each cluster also presented unique enriched pathways (Figure 6d-G and Figure S8).  Figure S8).
These results suggested that disc degeneration is not a uniform molecular process; but, similar to wide spectrum of degenerative morphological phenotypes , it is underscored by multiple molecular signatures.

| D ISCUSS I ON
While disc degeneration is complex and thought to have contributions from a multiple factors, aging is considered one of the major risk factors Ohnishi et al., 2018). Human disc degeneration has diverse phenotypes ranging from loss of disc height,  contribution of genetic background to these aging phenotypes is not well understood. Commonly used animal models with disc injury or genetic modifications have recapitulated some of these phenotypes, but do not provide insights why aging has such a wide spectrum of disc degeneration outcomes (Choi et al., 2018;Novais et al., 2019).
In this study, we show that the disc aging phenotype is dependent on the strain background in mice. While BL6 and SM/J mice present with a mild to fibrotic degeneration, respectively, LG/J mice, recapitulate several aspects of human age-dependent disc calcification and evidence aging transcriptomic changes that are both dependent and independent of the strain background.
BL6 mice are considered gold-standard model for aging studies and are widely used to characterize the aging of the spine Ohnishi et al., 2018). It is well known that human disc aging and degeneration comprises a wide spectrum of phenotypes, many of which are not recapitulated in the BL6 model, that is disc calcification (Chanchairujira et al., 2007;Cheung et al., 2009).
Interestingly, μCT analysis of LG/J mice showed ~80% prevalence of disc calcification in the caudal spine of old mice. Histological studies showed that this calcification was widespread and affected both the NP and AF. Additionally, presence of high concentration of free calcium not associated with a mineralization nodules and matrix proteins suggested its dystrophic nature. While studies of single gene knockouts affecting PPi metabolism showed ectopic mineralization of the AF and adjacent ligaments (Warraich et al., 2013;Zhang et al., 2016) LG/J mice recapitulated age-dependent disc calcification that affected both the NP and AF similar to elderly humans (Chanchairujira et al., 2007). Studies have also shown a causal relationship between systemic activity of TNAP that degrades PPi into Pi and ectopic calcification of soft tissues (Ziegler et al., 2017).
Interestingly, blood analysis showed that several mineral metabolism markers including free Ca 2+ and phosphorous were not affected by the strain background. While TNAP levels increased with aging in LG/J mice so did they in SM/J mice which also showed decreased levels of fetuin, a known inhibitor of ectopic calcification (Schäfer et al., 2003). Despite these differences in TNAP and fetuin, SM/J mice did not evidence disc mineralization suggesting that systemic mineral metabolism had little effect on the intradiscal mineralization.
Importantly, a previous study showed that LG/J mice despite having higher articular cartilage healing potential, evidenced increased susceptibility to synovial and meniscus calcifications following knee trauma (Rai, Schmidt, Hashimoto, Cheverud, & Sandell, 2015). Our results are in line with these findings that stressors like trauma and aging may promote mineralization in LG/J mice. Moreover, there was a higher incidence of calcification at Ca6-7, Ca7-8 and Ca8-9, whether increased mechanical/loading stress at these levels could be a contributing factor remains to be seen. Additionally, our results clearly showed that unlike caudal discs, lumbar discs of LG/J mice do not show calcification underscoring the importance of anatomic location, and consequently local mechanical environment, in predisposition to a degenerative phenotype (Teraguchi et al., 2014).
BL6, LG/J and SM/J discs showed different tissue and cellular organization with aging. Interestingly, BL6 caudal discs showed little increase in degeneration grades with aging as reported for lumbar and cervical spine (Ohnishi et al., 2018). This is not surprising as in humans, prevalence of disc degeneration is higher at lumbar L5/S1 and cervical C5/C6 than the other levels, suggesting that factors such as the local mechanical environment could play an important role in increasing the susceptibility to degeneration (Teraguchi et al., 2014). While SM/J mice presented higher grades of degeneration at 6M (Choi et al., 2018), they showed progressive changes with an increase in NP fibrosis. This was not unexpected since at 17 week, SM/J mice were shown to have fewer NP cells, and those remaining exhibited a hypertrophic chondrocyte-like phenotype responsible for secretion of collagenous fibrotic matrix (Choi et al., 2018). While most hypertrophic cells in NP did not survive at 23M, early changes in cell fate supports an overall decrease in aggrecan content (Choi et al., 2018). We also observed an accumulation of ARGXX neoepitope in the AF of old SM/J suggesting increased aggrecanase ac- and disc degeneration (Battié, 2009;Munir et al., 2018). Along with the higher grades of disc degeneration seen in SM/J at 6M, they differed the most in terms of gene expression compared to BL6 and LG/J. This suggested that severity and rate of progression of morphological changes were reflected in magnitude of changes at transcriptomic levels. Interestingly the genes and pathways that altered in SM/J such as cell differentiation, cell death, ion and proton transport and metabolic processes supported previous observations in this model (Choi et al., 2018). Surprisingly, although we found common DEGs between BL6 and LG/J with aging, most of them were divergent between the strains. Alterations in genes that are associ- osteoblasts. This idea is supported by a study that showed limited osteoblast-like characteristics by smooth muscle cells in the context of ectopic arterial calcification (Patel et al., 2019). Concerning increase in macrophage markers CD14, CD68 and CD163, studies have shown a strong link between macrophage activation with dystrophic calcification (Chen et al., 2016) and bone healing (Vi et al., 2015). Noteworthy, our recent results showed that NP cells possess immune cell-like properties and process extracellular antigens .
Histological analysis of human tissues have described various morphological features: disc fibrosis with acellularity, annular clefts, neovascularization, sclerosis of the subchondral bone; ectopic calcification of the endplate and/or the disc itself; herniation with increased senescent cells (Roberts, Evans, Trivedi, & Menage, 2006).
However, the contribution of these abnormalities in the etiopathogenesis of disc disorders as well as their actual cause-whether genetic predisposition, tissues response or altered mechanical environment-is not clear. To elucidate the genetic signatures of these phenotypes, we analysed the transcriptomic profiles of nondegenerated and degenerated human discs. Our analysis resulted in 7 distinct clusters of disc degeneration with unique molecular pathways. We found that several altered pathways in 6-month-old SM/J mice were represented in several human degenerated clusters: upregulation of cell death in 4 different clusters; cell differentiation, and response to stress in 3; and changes in metabolism and ion/intracellular transport in 3 clusters. In addition, inflammation and immune processes altered during BL6 aging were seen in some human clusters. Importantly, disc calcification is shown to be present in about ~6% of degenerated human discs (Chanchairujira et al., 2007). Interestingly, similar to enriched pathways in old LG/J mice, we found two of the human clusters showed enrichment of genes in pathways concerning response to stress, wound healing, cell death, endochondral bone, inflammation, cell division, phosphorous metabolism and extracellular matrix organization. These similarities support the hypothesis that common molecular pathways underscore disc degeneration in mice and humans. In conclusion, by analysing the aging NP signature in inbred mouse strains, as well as uncovering distinct transcriptomic profiles of human disc degeneration, our studies shed light on the wide spectrum of disorders presented in degenerative disc disease. Importantly, this study clearly demonstrated how genetic background can promote specific morphological phenotypes in the disc and highlights the utility of LG/J mice as a model to study intervertebral disc calcification.  (Staats, 1976), no systematic study has provided comparative estimates on cellular senescence activity between these three strains.

| Immunohistology and cell number measurements
De-paraffinized sections following antigen retrieval were blocked in 5% normal serum in PBS-T and incubated with antibodies against

| Microarray analysis and enriched pathways
Total RNA with RIN 5-9 was used for the analysis. Fragmented bi-

| Human microarray analysis
All sample details were reported by Kazezian et al. (2015) and available in GEO database, GSE70362 (Kazezian et al., 2015). Degenerative and nondegenerative clusters were obtained by hierarchical clustering of DEGs with a Euclidean distance < 0.5 cut-off. DEGs used for overrepresented biological processes enrichment analysis using PANTHER overrepresentation test, GO database annotation, binomial statistical test with FDR ≤ 0.05 of up-and downregulated DEGs between each degenerated cluster and cluster A (healthy).

| Micro-CT analysis
Micro

| FTIR imaging spectroscopy
IR spectral imaging data in the mid-IR region, 4,000-800/cm at 8/ cm spectral resolution and 25 μm spatial resolution was acquired from 10-μm-thick calcified caudal disc sections from 6M and 23M LG/J mice (n = 3 discs/group) using a Spectrum Spotlight 400 FT-IR Imaging system (Perkin Elmer, CT) (Choi et al., 2018;Mohanty et al., 2019). Absorbance for proteins in the amide I region, 1,665/cm; phosphate vibration region, 960 cm −1 and the carbonate at 870/cm were recorded (Berzina-Cimdina & Borodajenko, 2012). The preprocessed spectra were used for K-means cluster analysis to define anatomical regions and tissue types within the tissue section spectral images (Choi et al., 2018), which represent each analysed peak and ratio.
clustering images were obtained using Spectrum Image Software.

| Blood analysis
Blood was collected by heart puncture following sacrifice and centrifuged at 1,500 rcf, at 4 degrees, for 15 min to obtain plasma and send to measure analyte levels (IDEXX BioAnalytics). Fetuin-A was quantified using the mouse Fetuin-A/AHSG DuoSet ELISA (R&D Systems) according to the manufacturer's instructions.

| Statistical analysis
Analysis was performed using Prism7 (GraphPad, La Jolla). Data are represented as mean ± SEM. Data distribution was checked with Shapiro-Wilk normality test and the differences between two groups were analysed by t test or Mann-Whitney as appropriate.
The differences between 3 groups were analysed by ANOVA or Kruskal-Wallis for nonnormally distributed data followed by Dunn's multiple comparison test. Chi-squared test was used to analyse differences between distribution of percentages. p ≤ .05 was considered a statistically significant difference.

ACK N OWLED G M ENTS
This study is supported by the grants from the National Institute The authors would like to thank Drs. Andrzej Steplewski, William Querido and Nancy Pleshko for help with the FTIR.

CO N FLI C T O F I NTE R E S T
Nothing to disclosure.

DATA AVA I L A B I L I T Y S TAT E M E N T
Microarray data set that supports the findings of this study are openly available in GEO database, accession number GSE13 4955.