Abstract

   There are no appropriate mouse models to study the pathophysiology of
   spontaneous disc herniations in a wild-type setting. SM/J mice, a poor
   healer inbred strain, presented a high incidence of age-associated
   lumbar disc herniations with neurovascular innervations. Transcriptomic
   comparisons of the SM/J annulus fibrosus with human tissues showed
   shared pathways related to immune cell activation and inflammation.
   Notably, aged SM/J mice showed increased pain sensitization and
   neuroinflammation with altered extracellular matrix regulation in the
   dorsal root ganglia and spinal cord. There were increased T cells in
   the vertebral marrow, and cytometry by time-of-flight analysis showed
   increased splenic CD8+ T cells, nonspecific activation of CD8+ memory T
   cells, and enhanced interferon-γ production in the myeloid compartment.
   Single-cell RNA sequencing of peripheral blood mononuclear cells showed
   more B cells, with lower proportions of T cells, monocytes, and
   granulocytes. This study highlights the contribution of genetic
   background and aging to increased susceptibility of spontaneous
   intervertebral disc herniations in a clinically relevant murine model.
     __________________________________________________________________

   The SM/J mice highlight the role of genetics, aging, and neuroimmune
   axis in chronic back pain, a leading cause of disability.

INTRODUCTION

   With increased human lifespan, the interest in understanding
   age-dependent diseases has grown ([44]1). Musculoskeletal diseases,
   including osteoarthritis, osteoporosis, and intervertebral disc
   degeneration, are associated with aging, chronic pain, loss of
   mobility, and disability ([45]1). Disc degeneration is one of the major
   contributors to chronic low back pain (cLBP), a leading cause of years
   lived with disability ([46]2).

   While disc pathology manifests as a host of phenotypes, including
   fibrosis, calcification, and herniation ([47]3), understanding each is
   hampered by the lack of appropriate animal models. Although mouse
   models are developed to understand these degenerative phenotypes
   ([48]4–[49]6), studies of herniation are limited to injury models that
   capture local changes after acute herniations but fail to recapitulate
   age-related chronic changes ([50]7). Similarly, while transgenic and
   knockout mice have shown increased susceptibility to spontaneous
   herniations ([51]8, [52]9) and corresponding immune changes, they have
   limited translatability to humans because of single gene contributions
   and other comorbidities ([53]10). Notably, 85% of symptomatic human
   disc herniations are resorbed within 6 weeks because of
   neutrophil/macrophage and immune cell recruitment by the local
   pro-inflammatory and chemotactic signals ([54]11). Recent studies
   indicate that inhibiting these early inflammatory processes delays
   resorption and causes a transition to a chronic pain state ([55]12,
   [56]13).

   Many degenerative pathologies associated with chronic pain show
   increased innervation related to local expression of nerve growth
   factor (NGF) and protein gene product 9.5 ([57]14, [58]15). Notably,
   markers such as calcitonin gene-related peptide (CGRP) and isolectin B4
   (IB4) are used to discern peptidergic and nonpeptidergic nociceptive
   neurons ([59]16). While discs, facet joints, muscle, ligament, and
   fascia all can contribute to back pain, the sensory information is
   transmitted through identical circuits: local peripherical nociceptive
   neurons in dorsal root ganglia (DRG), the dorsal horn in the spinal
   cord, and lastly, the specific sensory brain region ([60]17).
   Furthermore, immune and glial activation in the central and peripheral
   nervous systems is thought to mediate the transition from acute to
   chronic pain ([61]18).

   In the present study, we provide comprehensive analyses of the
   pathophysiology and interdependence of spontaneous lumbar disc
   herniations and radiculopathy in SM/J mice, which we found are highly
   susceptible to age-dependent herniation pathology ([62]19). SM/J, a
   poor healer mouse strain, has become an attractive model for studying
   spontaneous disc degeneration because of the early onset of pathology
   in the caudal spine ([63]19, [64]20). These mice have a comparable
   lifespan to the commonly used C57BL/6 strain, and disc degeneration is
   not the result of a single gene mutation or external injury ([65]4). A
   similar degeneration process is also observed in genetically
   susceptible human populations, as suggested by twin studies ([66]21).
   Broadly, our studies highlight the strong association between immune
   system dysregulation in chronic pain because of disc herniation and
   emphasize the importance of genetic background to the susceptibility of
   age-associated herniations.

RESULTS

SM/J mice show a high incidence of age-associated lumbar disc herniations

   We characterized the lumbar spinal phenotype of SM/J mice at 6, 12, and
   20 months (M) and compared it to age-matched C57BL/6 mice. Histological
   analysis showed disc degeneration in SM/J at 6 M, which progressed in
   severity at 12 and 20 M, characterized by loss of nucleus pulposus
   (NP)/annulus fibrosus (AF) demarcation, decreased NP and AF
   cellularity, and increased AF clefts ([67]Fig. 1, A to C‴). While
   differences in grading distribution were smaller at 6 M, there was a
   clear increase in higher, more degenerative grades by 20 M ([68]Fig. 1,
   D and D′). As expected, both strains demonstrated age-associated lumbar
   disc degeneration; however, in SM/J, it was early onset and severe,
   evidenced by higher NP grades at 6 M and NP and AF grades at 12 and 20
   M ([69]Fig. 1, E and E′). Notably, ~90% of SM/J mice showed herniations
   at multiple lumbar levels by 20 M, with the highest incidence at L2-L5
   ([70]Fig. 1, F and G, and fig. S1A).

Fig. 1. SM/J mice show a high incidence of age-dependent disc herniation,
accompanied by cellular deficiencies and matrix remodeling.

   [71]Fig. 1.
   [72]Open in a new tab

   (A to C″) Histology of the 6-, 12-, and 20-M-old BL/6 and SM/J lumbar
   discs. SM/J shows altered tissue architecture and cell morphology,
   indicating disc degeneration in the NP, AF, and cartilage end plate
   (CEP). (D to E′) Modified Thompson grade distributions and averages
   show higher average scores in SM/J NP and AF tissues. (F and G)
   Prevalence of herniations. Lumbar levels from L1-S1: n[BL/6 6 M] = 5 (3
   females and 2 males); n[BL/6 12 M] = 6 (2 females and 4 males); n[BL/6
   20 M] = 9 (3 females and 6 males); n[SM/J 6 M] = 7 (4 females and 3
   males); n[SM/J 12 M] = 10 (6 females and 4 males); n[SM/J 20 M] = 8 (2
   females and 6 males). (H to J″) Staining and abundance of key markers
   of (H to H″) NP cell phenotype (CA3), (I to I″) local inflammation
   (IL-6), and (J to J″) cell viability (TUNEL). BL/6: 12 M (n = 5 to 6
   mice) and 20 M (n = 5 to 6 mice). SM/J: 12 M (n = 5 to 6 mice) and 20 M
   (n = 5 to 6 mice); two to three discs per animal. (K to M′) Picrosirius
   red staining and polarized light microscopy of lumbar disc sections. (N
   to N″) Analysis of percentage of thin (green), intermediate (yellow),
   and thick fibers (red). BL/6: 12 M (n = 5 to 6) and 20 M (n = 5 to 6).
   SM/J: 12 M (n = 5 to 6) and 20 M (n = 5 to 6); two to three discs per
   mouse. (O) Percentage of 20-M lumbar discs with collagen fibers in the
   NP. Kruskal-Wallis test for non-normally distributed data followed by
   Dunn’s multiple comparison test was used for (E) and (E′), presenting
   adjusted P value. A two-tailed t test or Mann-Whitney test was used for
   (H″) to (J″) and (N) to (N″), and chi-square test was used for (D),
   (D′), (G), and (O). Scale bars (5×), 200 μm. Scale bars (20×), 200 μm.

   We investigated the impact of herniations on cell phenotype, survival,
   and local inflammation. Unlike BL/6, SM/J mice did not show
   quantifiable staining of the NP phenotypic marker CA3 ([73]Fig. 1, H to
   H″) while showing an increased abundance of IL-6 in NP and AF ([74]Fig.
   1, I to I″). Furthermore, SM/J presented a higher percentage of
   terminal deoxynucleotidyl transferase–mediated deoxyuridine
   triphosphate nick end labeling (TUNEL)–positive cells, indicating
   compromised cell survival ([75]Fig. 1, J to J″). These results show
   that SM/J mice experience accelerated age-dependent disc degeneration,
   with increased incidence of herniations accompanied by local
   inflammation, cell death, and loss of the NP cell phenotype.

Aged SM/J mice evidence alterations in matrix composition and remodeling

   To understand extracellular matrix (ECM) health during aging and in the
   context of spontaneous herniations, collagen fiber thickness and
   markers of disc matrix integrity [collagen type II (COL2), cartilage
   oligomeric matrix protein (COMP), aggrecan proteoglycan (ACAN), and
   chondroitin sulfate (CS)] and degradation [ARGxx and matrix
   metallopeptidase 13 (MMP13)] were evaluated ([76]22). Picrosirius red
   staining showed a lower proportion of thin fibers and a higher
   proportion of intermediate and thick fibers in 6-M SM/J discs ([77]Fig.
   1, K, K′, and N to N″, and fig. S1B). While no differences were
   observed at 12 M, at 20 M, a higher proportion of thin fibers, as well
   as a lower proportion of intermediate fibers, with significantly higher
   NP fibrosis was noted ([78]Fig. 1, L′ to O, and fig. S1B). Notably, the
   collagen fiber thickness in SM/J discs remained relatively constant
   during aging; by contrast, in BL/6, a decrease in thin fibers and a
   concomitant increase in intermediate and thick fibers were noted,
   suggesting early-onset NP fibrosis in SM/J and age-associated AF
   collagen maturation in BL/6 ([79]Fig. 1, N and O).

   COL2 showed a marked reduction in abundance in NP and AF of 20-M SM/J
   mice (fig. S2, A to A″); however, COMP levels were comparable between
   the strains (fig. S2, B to B″). Notably, the abundance of ACAN, the
   primary disc proteoglycan, was lower in NP and AF compartments of SM/J
   mice (fig. S2, C to C″); however, the abundance of CS, a
   glycosaminoglycan side chain critical to the hydration of aggrecan and
   other sulfated proteoglycans, was decreased in the 20-M SM/J NP
   compartment but elevated in the AF at both 12 and 20 M, a
   characteristic feature associated with disc degeneration (fig. S2F).
   Moreover, a higher abundance of ARGxx, a maker of aggrecan degradation,
   and MMP13 was observed in SM/J discs (fig. S2, D to E″). These results
   indicate that SM/J herniations coincide with an altered ECM composition
   and increased catabolism.

AF tissue shows up-regulated DEGs associated with inflammatory activation
before herniation onset but diminished immune activity in the context of
persistent herniation

   To investigate the molecular signatures associated with increased
   susceptibility to disc herniation, microarray analysis of the AF tissue
   before and after the onset of herniations was performed using 12-M SM/J
   (degenerated, nonherniated discs) and BL/6 (healthy discs) and after
   herniation using 20-M SM/J (degenerated, herniated discs) and BL/6
   (degenerated, nonherniated discs) ([80]Fig. 2, A and B). We focused our
   analysis on lumbar AF tissue on the basis of our previous studies of
   caudal disc phenotype in SM/J mice (high NP degeneration; no
   herniation) and tumor necrosis factor–α (TNF-α) transgenic mice
   (healthy NP phenotype; higher susceptibility to disc herniation), in
   which the NP degeneration status did not predict the susceptibility to
   disc herniation ([81]4, [82]10). Principal components analysis showed
   distinct strain- and aging-based profiles ([83]Fig. 2B). The
   susceptibility of 12-M SM/J to developing hernia was associated with
   2595 differentially expressed genes (DEGs) [46% up-regulated, 54%
   down-regulated, false discovery rate (FDR) < 0.05, and fold change
   (FC) ≥ 2.0] (fig. S3C), wherein up-regulated DEGs were enriched for B
   and T cell activation, inflammation, apoptosis, p53, Janus
   kinase/signal transducers and activators of transcription signaling,
   and angiogenesis; and down-regulated DEGs showed enrichment for the
   regulation of metabolic processes, innate response, response to stress,
   and ion transport (fig. S3, D to F), highlighting possible contributors
   to AF disruption and eventual disc herniation. While both SM/J and BL/6
   up-regulated aging DEGs (FDR < 0.05 and FC ≥ 2.0) were enriched for T
   cell activation and apoptosis, SM/J aging down-regulated DEGs showed
   enrichment for the p53 pathway in contrast to BL/6 as well as
   cholecystokinin receptor and insulin/insulin-like growth factor (IGF)
   cascade (fig. S4A). To better understand the DEGs associated with the
   herniation phenotype independent of aging, we overlapped both up- and
   down-regulated DEGs of SM/J and BL/6 aging, uncovering unique genes in
   the SM/J group (herniation phenotype susceptibility) [602 (up) and 932
   (down)], common aging genes between the strains [2894 (up) and 2832
   (down)], and unique genes of the BL/6 group (fibrosis phenotype
   susceptibility) ([84]Fig. 2, D and E, and table S1). Herniation
   susceptibility DEGs were enriched for distinctive clusters of
   biological pathways in comparison to disc aging and fibrosis
   comparisons ([85]Fig. 2F and fig. S4, C and D). Up-regulated DEGs were
   associated with cell migration, positive regulation of
   phosphatidylinositol 3-kinase (PI3K) kinase, angiogenesis, collagen
   fibril organization, and shared axon guidance pathways with the
   fibrosis group. Down-regulated DEGs were associated with
   platelet-derived growth factor (PDGF) signaling, IGF cascade, PI3K, p53
   pathway, and interleukin (IL)/inflammation mediated by chemokine and
   cytokine pathways, the last two of which were also associated with
   up-regulated DEGs in the fibrosis groups. These results supported the
   hypothesis that herniation pathophysiology differs from normal aging
   and fibrosis phenotypes.

Fig. 2. Transcriptomic signatures of human AF degeneration share similarities
with SM/J mice.

   [86]Fig. 2.
   [87]Open in a new tab

   (A) Schematic summarizing study design for the transcriptomic analysis
   of AF tissues. (B) Transcriptomic profiles of 12-M BL/6 AF (n = 5
   mice), 20-M BL/6 AF (n = 4 mice), 12-M SM/J AF (n = 4 mice), and 20-M
   SM/J AF (n = 4 mice) tissues from mice clustered distinctly by
   principal components analysis. (C) Hierarchical clustering of SM/J
   aging DEGs (FC ≥ 2 and FDR ≤ 0.05). (D) Common and divergent
   up-regulated DEGs between SM/J AF and BL/6 AF aging. (E) Common and
   divergent down-regulated DEGs between SM/J AF and BL/6 AF aging. (F)
   Representative Panther Pathway processes from herniated group DEGs. (G)
   Human samples from [88]GSE70362 included in the comparative analysis
   with SM/J. (H) Hierarchical clustering of herniated versus degenerated
   DEGs (P ≤ 0.05). (I) Representative Panther Pathway processes from up-
   and down-regulated herniated versus degenerated DEGs. Panther analysis
   was performed using the PANTHER Overrepresentation Test with Panther
   Pathway and Panther Go Slim biological pathway database annotations and
   a binomial statistical test (FDR ≤ 0.05).

Transcriptomic signatures of human AF degeneration share similarities with
SM/J mice

   To evaluate the clinical relevance of SM/J transcriptomic changes, we
   compared them with a previously reported dataset ([89]GSE70362) of
   healthy, degenerated, and herniated human AF tissues ([90]Fig. 2G)
   ([91]23). The hierarchical clustering of these groups showed clear
   differences between healthy and degenerated/herniated tissues (fig. S5
   and S6). When degenerated and herniated samples were compared, a more
   complex relationship was revealed, with some degenerated samples more
   closely associating with herniated samples, reflecting that many
   herniated discs share features with discs exhibiting fibrotic
   degeneration ([92]Fig. 2H). Up-regulated DEGs (hernia versus
   degeneration group) showed an enrichment for PDGF, Ras, PI3K, FGF
   (fibroblast growth factor), EGF (epidermal growth factor) receptor, B
   cell, and IL-associated pathways ([93]Fig. 2I). By contrast,
   down-regulated DEGs were enriched for cell communication, glycolytic,
   antigen processing and presentation, sensory perception, and adaptive
   immune response ([94]Fig. 2I). Of note, SM/J unique DEGs also presented
   modulation of PDGF, PI3K, and IL and similar association of immune
   response with down-regulated DEGs, suggesting species-independent
   biological pathways correlated with the herniation phenotype. To
   further investigate this, we determined common DEGs between human
   (degeneration, no hernia) versus 12-M SM/J (degeneration, no hernia)
   and human (degeneration, with herniation) versus 20 M SM/J
   (degeneration, with hernia). Common DEG lists are presented in figs. S5
   to S7. Last, we used the CompBio tool to delineate commonality at the
   pathway level using the DEGs from the human herniated group and SM/J
   aging comparison (table S2) ([95]6). These results suggest a degree of
   conservation between human and SM/J AF degeneration, highlighting
   shared molecular pathways underpinning this pathology.

SM/J mice show altered structural integrity of subchondral bony end plate

   Sclerosis, increased bone mineral density, and decreased closed
   porosity of subchondral bone correlate with pain and local nerve growth
   in back pain and osteoarthritis models ([96]24). Because of this
   proposed relationship between the disc and the structural integrity of
   adjacent vertebrae and subchondral bony end plate, morphological
   metrics of the L3-S1 vertebral subchondral bony end plate were
   analyzed. 20-M SM/J mice showed an accelerated osteoporotic phenotype,
   evidenced by decreased cortical and trabecular thickness and number
   ([97]4). Subchondral bony end-plate volume did not differ between SM/J
   and BL/6 despite smaller vertebral size in SM/J ([98]Fig. 3, A to C).
   This finding coincided with higher tissue mineral density and cortical
   thickness of end plates in SM/J males and a higher bone volume fraction
   in females and males ([99]Fig. 3, D to F). In addition, SM/J males
   showed lower bone surface complexity and closed porosity ([100]Fig. 3,
   G and H).

Fig. 3. SM/J mice show an altered subchondral bone structure, tissue
innervation, and altered sensitivity to mechanical and temperature stimuli
with aging.

   [101]Fig. 3.
   [102]Open in a new tab

   (A and A′) Microcomputed tomography reconstruction of 20-M BL/6 and
   SM/J vertebrae, with the marked subchondral bone region of interest. (B
   to H) Analysis of subchondral bony end-plate parameters: (B) bone
   volume, (C) cross-sectional perimeter, (D) tissue mineral density, (E)
   cross-sectional thickness, (F) bone volume (BV/TV), (G) bone
   surface/volume (BS/BV), and (H) % closed porosity between 20-M BL/6 and
   SM/J males and females. n[BL/6 female] = 2; n[BL/6 male] = 3; n[SM/J
   female] = 2; n[SM/J male] = 2; four vertebrae (L3-L6 per mouse).
   Two-tailed t test or Mann-Whitney test was used as appropriate. (I to
   J″) Neurofilament and endomucin presence in BL/6 and SM/J discs at 12
   and 20 M, with the endomucin presence quantified. n = 6 to 8 mice per
   strain; two to three levels per mouse. Scale bars, 200 μm. (K to L″)
   Quantified NGF staining (K to L′) in 12- and 20-M BL/6 and SM/J NP (K″)
   and AF (L″) tissues. (M to O) Sensitization analysis: (M) mechanical
   (Von Frey), (N) hot (Hargreaves), and (O) cold stimulus. n[BL/6 12
   M] = 6 to 11; n[BL/6 20 M] = 8; n[SM/J 12 M] = 11 to 14; n[SM/J 20
   M] = 6 to 7. (P) Schematic of gait analysis and parameters. (Q to X)
   Gait parameters: (Q) swing time, (R) paw angle variability, (S) stance
   width variability, (T) step angle, (U) stride length variability, (V)
   swing duration, (W) paw area peak stance, and (X) gait symmetry. n[BL/6
   12 M] = 9; n[BL/6 20 M] = 6; n[SM/J 12 M] = 8; n[SM/J 20 M] = 7.
   Kruskal-Wallis test for non-normally distributed data followed by
   Dunn’s multiple comparison test was used for (M) to (O), with adjusted
   P value. Two-tailed t test or Mann-Whitney test was used to compare
   differences between two groups as appropriate for (B) to (H), (K′) to
   (L″), and (K) to (X). Chi-square test was used for (J″). (I to J′)
   Scale bars (20×), 200 μm. (K to L′) Scale bars (10×), 200 μm.

SM/J mice show peripheral vascularization and innervation of the herniated
discs

   Disc herniation is tightly correlated with cLBP and radiculopathy.
   Neurofilament, a marker of axonal ingrowth ([103]16), and endomucin, a
   marker of neovascularization commonly associated with neoinnervation
   ([104]15), were used to investigate nerve and blood vessel ingrowth in
   herniated SM/J discs. We observed that endomucin- and NF-stained
   structures were in proximity outside the disc at 12 M but were present
   within the outer AF of herniated discs in 20-M SM/J ([105]Fig. 3, I to
   J″). Endomucin staining was pronounced in more 20-M-old SM/J mice
   compared to age-matched BL/6 controls, suggesting increased
   vascularization in this group ([106]Fig. 3J). These results confirm an
   association between disc herniation and local axonal and blood vessel
   ingrowth, as reported in humans ([107]25). In addition, neuronal and
   pain marker NGF was significantly higher in 20-M SM/J NP and AF tissues
   than in 20-M BL/6 ([108]Fig. 3, K to L″), suggesting that herniations
   in SM/J mice may stimulate neuropathic pain.

SM/J mice show increased sensitization to mechanical and cold stimuli, along
with changes in gait

   Human herniation is often associated with several neuropathic symptoms,
   including paresthesia, allodynia, hyperalgesia, and heat sensitization
   ([109]26, [110]27). To assess similar behavioral changes, mechanical
   (Von Frey filament) and temperature (Hargreaves thermal; cold plantar
   assay) sensitization was measured in SM/J and BL/6 mice ([111]28,
   [112]29). Compared to BL/6, SM/J mice demonstrated a lower threshold
   for mechanical stimulation inducing a pain-like response during aging
   ([113]Fig. 3M). In addition, SM/J showed a decreased latency for a
   noxious cold stimulus inducing a pain-like response than BL/6 at 12 and
   20 M ([114]Fig. 3N). By contrast, Hargreaves testing (heat sensitivity)
   showed increased desensitization in 12-M SM/J mice, which persisted at
   20 M ([115]Fig. 3O). To study the correlation of gait features with
   neurological pathologies or locomotion-related pain, gait analysis was
   performed ([116]Fig. 3P) ([117]30). SM/J mice presented an
   age-dependent increase in swing time; however, compared to BL/6, the
   swing time was shorter in SM/J, independent of age ([118]Fig. 3Q). Paw
   angle variability increased with age in BL/6, and a similar trend was
   observed in SM/J mice ([119]Fig. 3R). While the stance width
   variability rose with age in BL/6, it remained lower than in SM/J at 12
   M ([120]Fig. 3S). There were no significant changes in other parameters
   associated with ataxia (step angle and stride length variability),
   impaired mobility (swing duration), limb injury (paw area peak stance),
   or gait symmetry ([121]Fig. 3, T to X). These findings show that SM/J
   replicates clinical behaviors commonly seen in human patients with disc
   herniation and radiculopathy.

SM/J DRGs show enrichment of DEGs associated with the glial cells,
neuroinflammatory response, cell-matrix adhesion modulation, and immune
system regulation

   To study the molecular changes in DRGs concordant with herniation and
   behavioral changes, transcriptomic analysis of lumbar DRGs was
   performed. While SM/J profiles diverged with age, BL/6 profiles did not
   mirror this separation ([122]Fig. 4A). Source-of-variation analysis
   showed that the compounded effects of age and strain (“Condition”) were
   the highest source of variation among samples, whereas sex-based
   effects were negligible in both strains ([123]Fig. 4B). Hierarchical
   clustering of 12-M SM/J and BL/6 showed divergent transcriptomes
   between the strains (fig. S8A). There were 188 up-regulated DEGs
   (FDR < 0.05 and FC ≥ 2), which did not enrich in any biological
   pathways (PantherDB) (fig. S8, B and C), whereas 167 down-regulated
   DEGs were associated with the response to interferon-β (IFN-β),
   positive regulation of response to a biotic stimulus, activation of the
   innate immune response, antigen processing, and presentation of peptide
   antigen pathways (fig. S8, B and C). Comparison of DEGs from 20-M SM/J
   with 12-M SM/J showed distinct clustering during aging (fig. S9A), and
   a total of 401 DEGs was noted (327 up-regulated and 74 down-regulated)
   (fig. S9B). There was no biological enrichment within the
   down-regulated DEGs; however, up-regulated DEGs were associated with
   ECM organization and adhesion regulation; modulation of IL-1, IL-4,
   IL-6, IL-10, and IL-12; regulation of IFN-α and transforming growth
   factor–β pathways; synapse pruning; cellular response to amyloid-β;
   regulation of neuron death; astrocyte activation and microglial cell
   migration; Toll-like receptor; and negative regulation of B cell
   receptor signaling (fig. S9, B and C).

Fig. 4. DRGs in SM/J mice show the enrichment of DEGs associated with the
glial cells, neuroinflammatory response, cell-matrix adhesion modulation, and
immune system regulation.

   [124]Fig. 4.
   [125]Open in a new tab

   (A) Principal components analysis of transcriptomic profiles of 12- and
   20-M BL/6 (n = 3 females and 3 males per time point) and SM/J (n = 3
   females and 3 males per time point) DRG tissues, (B) weighted source of
   variation, and (C) volcano plot showing up- and down-regulated DEGs
   from the 20-M SM/J DRGs versus 20-M BL/6 DRG comparison (FC ≥ 2 and
   FDR ≤ 0.05). (D) Pie chart showing the distribution of up- and
   down-regulated DEGs and enriched GO processes from the 20-M SM/J versus
   BL/6 DRGs. (E) Venn diagram of common up-regulated DEGs in 20-M SM/J
   versus 12-M SM/J and 20-M SM/J versus 20-M BL/6 but not present in 12-M
   SM/J versus 12-M BL/6. (F) Pathway enrichment analysis of 115 DEGs
   identified in (E); analysis performed in PantherDB using the
   statistical overrepresentation test with Gene Ontology database
   annotations and a binomial statistical test (FDR < 0.05). Values in
   parentheses demonstrate the number of genes within the annotation set
   relative to the predicted number of genes. (G to J″″) Quantitative
   immunohistochemistry of the dorsal horn of the spinal cord in the
   lumbar region: (G to G″″) IB4, (H to H″″) CGRP, (I to I″″) GFAP, and (J
   to J″″) IBA1. A.U., arbitrary units. n[BL/6 12 M] = 6; n[BL/6 20
   M] = 4; n[SM/J 12 M] = 6; n[SM/J 20 M] = 4; two spine levels per
   animal. Two-tailed t test or Mann-Whitney test was used as appropriate.
   Nonparametric Spearman correlation test was used for analysis in (D) to
   (G). Scale bars (A to H′), 200 μm.

   Hierarchical clustering of 20-M SM/J versus BL/6 showed clear
   strain-based segregation of DEGs (fig. S9A) with 364 up-regulated and
   431 down-regulated DEGs in SM/J ([126]Fig. 4D). Up-regulated DEGs (fig.
   S10B) showed enrichment of negative regulation of cell-matrix adhesion
   and hyaluronan metabolic process; regulation of IL-1β, IL-6, and IL-8;
   transforming growth factor–β activation; response to axon injury;
   cellular response to serotonin; glial cell migration; sensory
   perception of pain signaling; regulation of apoptotic cell clearance;
   ERK1 (extracellular signal–regulated kinase 1) and ERK2 cascade;
   modulation of Toll-like receptor signaling; macrophage inflammatory
   protein 1α production; and negative regulation of T cell proliferation
   ([127]Fig. 4D). Down-regulated DEGs (fig. S10B) were enriched for
   trans-synaptic signaling, sensory perception of pain, and other
   signaling pathways ([128]Fig. 4D).

   To identify any DEGs associated with the pain-related behaviors
   exhibited by 20-M SM/J and to minimize the effect of age- and
   strain-based variation, we investigated the DEGs common between the
   SM/J aging and 20-M SM/J versus BL/6 comparisons and subtracted the
   baseline strain-associated DEGs identified from 12-M SM/J versus BL/6
   comparison ([129]Fig. 4E and fig. S1L). Eight down-regulated
   transcripts, Gp1bb, Spry3, Izumo4, Ephx4, Plpp4, Foxj1, Prrt2, and
   Arhgdig (fig. S11, A and B), and 115 up-regulated DEGs met these
   criteria ([130]Fig. 4E). Consistent with the observed pathophysiology
   of SM/J mice in this study, these DEGs were associated with the ECM,
   cytokines, nervous system, cell viability, and immune regulation
   pathways ([131]Fig. 4F). These results provide new insight into the
   pathological processes in DRGs in this novel model of spontaneous disc
   herniation.

SM/J mice show higher staining for IB4 and GFAP in the lumbar dorsal horn

   We evaluated the abundance of IB4-binding (nonpeptidergic) and
   CGRP-positive (peptidergic) nociceptive axons in the dorsal horn of the
   thoracic (adjacent to nonherniated discs) and lumbar (herniated) spinal
   cord ([132]16). In the lumbar region, an abundance of IB4-positive
   axons was higher in SM/J mice at 12 and 20 M, whereas the abundance of
   CGRP-positive axons was maintained across all groups ([133]Fig. 4, G to
   H″″). Microglia and astrocytes play essential roles in pain modulation,
   and the reactivity of these cell types in the spinal cord dorsal horn
   is associated with multiple chronic pain states ([134]17, [135]31). No
   differences were observed in levels of the microglial marker Iba1, but
   20-M SM/J showed higher levels of the reactive astrocytic marker glial
   fibrillary acidic protein (GFAP) ([136]Fig. 4, I to J″). This increase
   in GFAP abundance coincided with a broader spatial distribution of
   astrocytes in the deeper layers of the dorsal horn, something not
   observed in BL/6 or 12-M SM/J mice. By contrast, no significant changes
   were observed in these markers in the thoracic region (fig. S11, D and
   E). These findings suggest a possible role of IB4 nociceptive neurons
   and astrocytes in pain sensitization in SM/J mice with herniations.

SM/J mice show altered plasma cytokine levels and increased T cell and
neutrophil recruitment in the vertebral marrow

   We assessed the impact of strain and aging on the systemic burden of
   inflammation by evaluating the abundance of pro- and anti-inflammatory
   markers. Plasma levels of IL-10, IL-15, IL-17A, IL-18A/F, IL-1β,
   IL-27p28/IL-30, monocyte chemoattractant protein 1, and macrophage
   inflammatory protein–3α (MIP-3α) rose during aging in BL/6 ([137]Fig.
   5, A to R, and fig. S11). While SM/J showed a similar increase in
   IL-17A and MIP-3α levels with aging, an increase in IL-16, IL-33, and
   MIP-1α was also noted ([138]Fig. 5, P and Q, and fig. S12). 12-M SM/J
   had higher baseline levels of IL-21, IL-23, and IL-6 than BL/6
   ([139]Fig. 5, H, I, and L, and fig. S12). At 20 M, SM/J evidenced lower
   IL-2, IL-4, and IFN-γ compared to BL/6. Levels of IL-17C, IL-22, IL-5,
   IL-9, IP-10 (IFN-γ–induced protein 10), keratinocyte
   chemoattractant/human growth-regulated oncogene (KC/GRO), and MIP2 were
   consistent among all groups. Notably, IFN-γ, IL-1β, IL-2, and IL-4
   showed a negative correlation to sensitization to mechanical stimuli
   ([140]Fig. 5, S to V). We evaluated the recruitment of B cells (CD19),
   T cells (CD3), neutrophils (Ly6G), macrophages (F4/80), and activated T
   cells and macrophages [major histocompatibility complex class II
   (MHCII)] in the vertebral marrow ([141]32). While no changes were
   observed in CD19, F4/80, or F4/80-MHCII double-positive cells, CD3 and
   MHCII abundance was higher in 20-M SM/J mice than in age-matched BL/6
   and 12-M SM/J mice ([142]Fig. 5, W to Z). Higher levels of Ly6G in 20-M
   SM/J mice relative to 20 M BL/6 mice were also observed ([143]Fig. 5, W
   to Z). These studies support the idea that the immune system
   contributes to disc herniation–associated radiculopathy.

Fig. 5. SM/J mice show altered plasma cytokine levels and increased T cell
and neutrophil recruitment in the vertebral marrow.

   [144]Fig. 5.
   [145]Open in a new tab

   (A to R) Multiplex analysis of pro-inflammatory, dual function, and
   anti-inflammatory proteins in plasma from 12- and 20-M BL/6 and SM/J
   mice. t test or Mann-Whitney test was used as appropriate. (S to V)
   Pearson correlation between INF-γ, IL-1, IL-2, and IL-4 plasma levels
   and Von Frey sensitivities. n = 4 to 6 mice per strain per time point.
   (W to Z″″″) Quantitative immunohistochemistry of (W to W″″) CD19, (X to
   X″″) CD3, (Y to Y″″) Ly6G, and (Z to Z″″) MHCII and F4/80 (individually
   and ratio quantification) at 12 and 20 M in the vertebral bone (n = 6
   to 10 vertebrae per strain per time point). Two-tailed t test or
   Mann-Whitney test was used as appropriate. Scale bars (A to H′), 100
   μm.

SM/J mice show altered immune cell profiles in peripheral circulation

   Immune cells are key modulators of cLBP development ([146]13). To
   investigate the circulating immune cell profile and understand their
   association with pain pathology in SM/J mice, we performed a
   single-cell RNA sequencing (scRNA-seq) of peripheral blood mononuclear
   cells (PBMCs) ([147]Fig. 6A). After normalization, the cells were
   grouped into 21 clusters according to the FindNeighbors and
   FindClusters functions within the Seurat single-cell sequencing
   analysis package ([148]Fig. 6A). Furthermore, the data were visualized
   according to age and strain and cell types were labeled using the
   SingleR annotation tool, demonstrating ([149]33) marked changes in cell
   populations across strains and with aging ([150]Fig. 6, B to C″, and
   fig. S9A). No major changes in the distributions of cell types were
   observed across strains at 12 M, but by 20 M, SM/J showed lower
   proportions of monocytes and T cells and higher proportions of B cells
   and macrophages than BL/6 ([151]Fig. 6, C to C″).

Fig. 6. SM/J mice show altered immune cell profiles in peripheral
circulation.

   [152]Fig. 6.
   [153]Open in a new tab

   (A) Schematic showing experimental design and unbiased clustering of
   PBMC transcripts. n[12 M BL/6] = 2 females and 2 males; n[12 M
   SM/J] = 2 females and 1 male; n[20 M BL/6] = 2 females and 2 males;
   n[20MSM/J] = 4 males. (B) Clustering of PBMC transcripts according to
   strain and time point. (C to C″) Proportion of cell populations labeled
   with the SingleR program. (D and E) UMAP labeled according to mouse
   strain and age, showing cell populations identified to be unique to
   SM/J at 20 M, when mice show disc herniations and pain behaviors: (D) B
   cells and (E) macrophages. (D′) GSEA in B cells based on DEGs
   identified in 20-M SM/J PBMCs. (D″) Violin plots of most DEGs in B
   cells in the SM/J pain group. (E′) GSEA of macrophages from 20-M SM/J.
   (E″) Violin plots of most DEGs in macrophages in the SM/J pain group
   (n ~ 5000 cells per strain per time point). Chi-square test was used
   for (C″).

   Previous studies have suggested the dependence of immune cell
   heterogeneity on inflammatory versus quiescence state, function, and
   origin ([154]34, [155]35). To visualize heterogeneity, cell populations
   were subset and visualized with uniform manifold approximation and
   projection (UMAP), showing distinct differences in cell populations
   according to pain status–specific differences—B cells and macrophages
   ([156]Fig. 6, D and E, and fig. S13, B and C)—and strain—T cells
   (subdivided into CD4+ and CD8a+ T cells), monocytes, natural killer
   (NK) cells, and granulocytes (fig. S14). Gene set enrichment analysis
   (GSEA) was conducted on DEGs defined by P < 0.05 for each cell
   population using the clusterProfiler package ([157]36). In B cells of
   20-M SM/J, GSEA showed the activation pathways related to T cell
   tolerance induction, Ras protein signal transduction, and immune system
   processes ([158]Fig. 6F′ and fig. S13, B′ and B″). On the other hand, B
   cells showed suppression of negative regulation of regulatory T cell
   differentiation and purine biosynthesis process, and top DEGs included
   Bach2, Cmss1, Ripor2, Ebf1, and Cdk8 ([159]Fig. 6, D′ and D″, and fig.
   S14, B′ and B″). While the proportion of macrophages in the pain group
   was higher, their DEGs showed suppression of defense- and immune
   response–related pathways and activation of platelets, blood
   coagulation, and serotonin uptake, with top DEGs including Nrgn,
   Clec1b, Itga2b, Prg4, and Chil1 ([160]Fig. 6, E′ and E″, and fig. S14,
   E′ and E″). SM/J CD4+ T cells showed DEGs associated with intracellular
   molecular regulation–related pathways (fig. S14B′). Despite being fewer
   in the 20-M SM/J group, CD8+ T cells were enriched for the activation
   of an innate response, positive regulation of leukocyte-mediated
   cytotoxicity and immunity, and cellular response to IFN-γ and for
   suppression of negative regulation of Wnt signaling, translation, and
   complement activation (fig. S14C′). GSEA of SM/J monocytes showed
   activation of positive regulation of T cell activation, antigen
   processing and presentation, and lymphocyte-mediated immunity and
   suppression of response to a stimulus and regulation of biological
   processes (fig. S13E′). Last, SM/J NK cells were enriched for the
   activation of the innate immune response, NK cell–mediated immunity,
   and apoptotic cell clearance (fig. S14F′).

   To explore the agreement between scRNA-seq and plasma analyses, we
   looked for cell populations expressing Ifng. We found T cells and NK
   cells to be the two predominantly Ifng-expressing populations, and
   trends in expression, considered with the cell population trends,
   confirm the reduced overall expression in 20-M SM/J mice (fig. S15, A
   to A″). These analyses suggest that dysregulation of the immune system,
   characterized by suppression of basic immune and regulatory processes
   in multiple cell types, of SM/J mice contributes to chronic disc
   herniation and pain.

SM/J mice show an altered immune cell profile in the spleen

   To further interrogate the associations between discogenic and
   radicular pain and the systemic immune system modulation, cytometry by
   time-of-flight (CyTOF) analysis was conducted on splenocytes from 12-
   and 20-M mice using 16 cell-surface markers and 8 intracellular
   cytokines (20 M only); a previous study has shown that the splenocyte
   profile reliably correlates with the systemic adaptive immune response
   ([161]37). Using established immune cell signature profiles, 15 cell
   populations were identified on the basis of an initial gating for CD45+
   cells ([162]Fig. 7A). This cell identification strategy also coincided
   with our scRNA-seq data, aligning with several populations identified
   in that analysis (fig. S16B). At 12 M, SM/J showed a significantly
   higher frequency of T cells and NK cells than BL/6 ([163]Fig. 7B and
   fig. S16, A to A″″), whereas at 20 M, this trend reversed for NK cells
   and there was no difference between two strains for overall T cell
   numbers. However, among T cells, there was a higher frequency of the
   CD8+ T cells in SM/J splenocytes ([164]Fig. 7, B′ to C′, and figs. S16,
   B to B″″, and S11). During aging, the NK cell frequency rose, and CD8+
   T cell and myeloid cell frequencies decreased in BL/6 ([165]Fig. 7C and
   fig. S16, C to C″″). Frequencies of NK cells, CD4+ memory T cells, and
   CD4+ naïve T cells decreased during aging in SM/J ([166]Fig. 7, C′ to
   I, and fig. S16, D to D″″).

Fig. 7. CyTOF analysis of splenocytes reveals age- and strain-based frequency
differences in immune cell populations, and intracellular marker abundance in
20-M SM/J splenocytes indicates possible repression of the innate immune
system and overactivation of CD8+ T cell populations.

   [167]Fig. 7.
   [168]Open in a new tab

   (A) Schematic showing the workflow for CyTOF analysis of splenocytes
   from BL/6 and SM/J mice and the cell population labeling hierarchy.
   n[12 M BL/6] = 2 females and 2 males; n[12 M SM/J] = 2 females and 1
   male; n[20 M BL/6] = 2 females and 2 males; n[20 M SM/J] = 4 males. DC,
   dendritic cell. (B to C′) Volcano plots showing cell populations in
   strain- and age-based comparisons based on FDR and FC. (D to I)
   Frequency plots showing age- or strain-based differences in (D) NK
   cell, (E) T cell, (F) CD8+ T cell, (G) myeloid cell, (H) CD4+ memory T
   cell, and (I) CD4+ naïve T cell populations. (J) Heatmap of splenocytes
   from 20-M SM/J versus BL/6 showing intracellular cytokine and
   pro-inflammatory marker abundance within each cell population
   identified through CyTOF analysis. (K to Q) Significant changes in
   cytokine/pro-inflammatory marker abundance in (K) myeloid cells, (L)
   conventional dendritic cells, (M) plasmacytoid dendritic cells, (N)
   Lin− cells, (O) CD8+ naïve T cells, (P) CD8+ T cells, and (Q) CD8+
   memory T cells. Two-tailed t test or Mann-Whitney test was used as
   appropriate.

   To investigate possible functional differences in the immune cell
   populations in mice with herniated discs, intracellular cytokine levels
   were determined at 20 M ([169]Fig. 7J and fig. S17). At 20 M, SM/J mice
   showed significantly higher numbers of myeloid cells producing IFN-γ
   ([170]Fig. 7K). These inflammatory myeloid cells that include
   granulocytes, monocytes, and macrophages (M1 macrophages) depicted that
   an inflammatory type 1 immune response [T helper 1 (T[H]1) response]
   characterized by IFN-γ production was predominant in SM/J mice. Lower
   production of anti-inflammatory T[H]2 cytokines IL-4 and IL-5 in
   plasmacytoid dendritic cells ([171]Fig. 7L) and CD8+ naïve T cells
   ([172]Fig. 7M), respectively, corroborated the inflammatory
   environment. IL-2 was also lower in conventional dendritic cells
   ([173]Fig. 7N). Furthermore, nonspecific activation of CD8+ T cells
   ([174]Fig. 7O) and CD8+ memory T cells ([175]Fig. 7P) was also
   observed. Overall, there was enhancement of an inflammatory T[H]1
   response and suppression of an anti-inflammatory T[H]2 immune response.
   The T[H]17 immune response was also suppressed, as seen by
   significantly lower production of IL-17A by Lin− cells ([176]Fig. 7Q)
   as well as myeloid cells ([177]Fig. 8A), indicating that an
   inflammatory T[H]1 immune response, not T[H]17, underscored the
   phenotype. Significant changes in cytokine profiles of other cell
   populations were not observed (fig. S12).

Fig. 8. Summary schematic describing the SM/J herniation phenotype.

   [178]Fig. 8.
   [179]Open in a new tab

   Working hypothesis and phenotypic changes found in the NP, AF,
   subchondral bone, DRGs, spinal cord, plasma, and immune tissues of SM/J
   mice relevant to herniation and chronic pain–related behaviors.

DISCUSSION

   Intervertebral disc degeneration comprises multiple degenerative
   phenotypes ([180]3). Recently, SM/J mice were shown to evidence
   early-onset, spontaneous disc degeneration characterized by advanced
   fibrosis in the caudal spine during aging compared to BL/6 and LG/J
   strains that show milder fibrosis and dystrophic calcification,
   respectively ([181]4, [182]19). These studies highlight the essential
   contribution of genetic background to degenerative processes,
   reinforcing findings from human twin studies that show that familial
   background is a larger contributor to disc degeneration than aging and
   mechanical loading ([183]21). This study shows that SM/J mice are
   highly susceptible to age-associated lumbar disc herniations. These
   results suggest that while genetics govern susceptibility to disc
   degeneration, anatomical/mechanical factors may influence and drive the
   final degenerative outcome ([184]4). This idea aligns with studies of
   LG/J mice that show age-dependent calcification restricted to the
   caudal spine ([185]4) and a significantly higher prevalence of disc
   degeneration at L3-5 of the human spine than the other levels ([186]3).
   However, the effect of genetic background and age as independent
   variables in this study must not be overlooked. The SM/J susceptibility
   to disc herniation aligns with previous genome-wide association studies
   in this model and the Twin Spine Study, suggesting genetic background
   as the primary modulator of regeneration potential and disc
   degeneration process ([187]21, [188]38). In humans, it is acknowledged
   that spinal lordosis and unit rigidity make L3-4 and L4-5 levels more
   susceptible to herniation, but with aging and a consequent decrease in
   lordotic angle, the prevalence increases in the upper lumbar levels
   ([189]39, [190]40). In line with these clinical findings, our study
   describes a model of spontaneous disc herniation associated with aging,
   which may explain a higher herniation frequency in upper lumber levels.
   While we cannot state with certainty the precise contribution of
   mechanics or lordosis in SM/J mice, it is important to highlight that
   the mechanics, tissue structure, and even embryonic signaling change
   across the mouse spine, which precede the regional susceptibilities
   across the spine ([191]41–[192]43). It is also true that multilevel
   disc herniation is less prevalent in humans, with sources reporting 53
   to 80% of herniations as a single level, 12 to 22% as two levels, and
   fewer presenting more than three affected levels ([193]44, [194]45).
   However, these generalizations may include a selection bias, as studies
   evaluating the thoracic region and older patients show increased
   incidence of multilevel disc herniation ([195]46, [196]47).
   Unfortunately, there are very little data on the prevalence and
   distribution of intervertebral disc herniation in the elderly
   population. In this sense and understanding of the impossibility of an
   animal model fully capturing every feature of human disease, SM/J mice
   present a unique opportunity for study, recapitulating the increased
   susceptibility to disc herniation with aging. Of note, 20-M-old mice
   correspond to elderly humans aged about 60 to 65 years old. In
   addition, previous work by Rajasekaran et al. ([197]48) suggested that
   disc herniation occurs more from end-plate junction failure than AF
   rupture. While there is a limited ability to assess AF rupture versus
   end-plate failure in SM/J, and SM/J mice may be a more severe
   herniation case than the average human patient, this model successfully
   mimics the human phenotype in the sense that SM/J herniation is not
   induced by acute injury or single-gene manipulation, as has
   historically been the case in mouse models of disc herniation.

   Considering that the ECM plays a critical role in disc function, the
   lack of age-associated thickening of collagen fibers in SM/J suggested
   diminished remodeling of the AF, possibly due to low cell density
   ([198]49). On the other hand, an increased abundance of MMP13 and ARGxx
   and decreased ACAN with increased CS levels suggested the breakdown of
   the aggrecan-rich matrix in SM/J mice, with possible compensatory
   efforts by other sulfated proteoglycans, such as versican ([199]50).
   This loss of ECM integrity likely causes loss of disc structural
   integrity and mechanical properties ([200]19). Notably, ACAN levels and
   turnover directly modulate axonal ingrowth and pain in degenerated
   human discs and arthritis models ([201]51, [202]52). In addition, the
   increased CS levels in the inner AF of SM/J mice may result in a loss
   of demarcation between the NP and AF and increased transition zone
   instability, further exacerbated by reduced abundance of COL2,
   promoting AF weakness and consequently increasing susceptibility to
   herniation ([203]53). Furthermore, reduced ACAN and CS levels with
   aging are likely an important contributor to a higher propensity for
   peripheral nerve ingrowth in herniated SM/J discs ([204]54, [205]55).
   When AF DEGs from 12-M SM/J were cross referenced with those from
   humans susceptible to disc herniation, we observed a down-regulation of
   ECM-related transcripts, further supporting a key role of ECM
   homeostasis in preventing herniations and previous findings showing
   fewer, thinner, and more disorganized collagen I fibers in human
   herniated discs ([206]8, [207]10, [208]56).

   Comparing SM/J findings to those of degenerated human AF ([209]57,
   [210]58) showed that dysregulation in cell cycle, cell viability,
   extracellular structure, and neutrophil activation pathways was shared
   between the two species. Previous studies have indicated the importance
   of the immune system in tissue healing after herniation ([211]59), and
   while the disc compartment is immune privileged ([212]49), our results
   suggest that dysregulation of local innate immune response may precede
   acute herniation and prevent efficient resorption of the herniated mass
   and scar tissue formation. Similarly, human TNF-α transgenic mice
   showed higher susceptibility to disc herniation associated with
   compromised AF integrity and local alterations in immune cell
   recruitment ([213]10, [214]60). Notably, despite higher vertebral
   osteopenia in old SM/J ([215]4), these mice show higher subchondral
   bone thickening. These results align with earlier observations, wherein
   an aberrant increase in subchondral bone volume/thickness in lumbar
   vertebrae and knee joints promoted by osteoclastic secretion of
   netrin-1 induces sensory innervation and pain ([216]14, [217]24).

   cLBP is the primary clinical symptom associated with disc herniations
   ([218]61). In addition, patients often experience radiating pain in
   their limbs, increased sensitization to mechanical and thermal stimuli,
   and altered gait ([219]61)—phenotypes we noted in aged SM/J mice.
   Before our investigation, the only behavioral study of SM/J mice
   evaluated depression-like behavior in 4- and 8-week-old animals and
   showed a lower response during tail suspension and forced swim tests
   ([220]62, [221]63). While our studies did not investigate depression as
   a possible contributor to pain, it is critical to note that the prior
   studies on depression in SM/J mice evaluated juvenile mice, as opposed
   to aged mice, and in the present study, SM/J mice with pain behaviors
   did not display other indicators of depression, such as fasting, weight
   loss, or altered grooming ([222]4, [223]64). Some studies in the
   literature suggest a sex-dependent influence on disc degeneration and
   back pain ([224]65). However, in our current study, we did not observe
   a significant impact of sex on the degenerative phenotypes. This may be
   attributed to the stronger influence of genetic background, which could
   overshadow any potential sex-related effects.

   Notably, while we cannot completely exclude a possible contribution of
   degeneration of other joints, namely ankle, knee, and hip, in our
   behavioral findings, these pathologies usually display a longer swing
   phase, slower swing speed, and decreased paw area in gait analysis,
   which we did not observe in our cohorts. Together with the changes
   observed in the lumbar region, but not at the thoracic region of the
   spinal cord, which corresponded to the lumbar levels of disc hernias
   and behavioral changes (leg innervation region), our data suggest a
   higher contribution of the spine phenotype to the neuropathic behavior.
   As in humans, the ideal method to characterize disc herniation is by
   magnetic resonance imaging studies; while this was not feasible, by
   performing sequential transverse and coronal sections, we were able to
   capture the herniated tissue through AF, speculating that in the most
   severe case, this herniated material may ingress into the foramen or
   the spinal canal, promoting local inflation, culminating in the
   behavioral changes we observed in our study.

   Mechanistically, degenerated SM/J discs had higher levels of IL-6, also
   noted in disc puncture models and osteoarthritis ([225]66, [226]67),
   which serves as a chemoattractant for axonal ingression, promoting
   local neuropathic pain ([227]68). Accordingly, increased vascular and
   sensory innervations were observed in the outer AF but, unlike other
   studies, not in the NP of herniated SM/J discs ([228]69). This finding
   may be explained by hernia remodeling stages: acute versus absorbed
   versus chronic ([229]5, [230]10). In humans, herniated discs show the
   presence of CGRP- and IB4-positive nociceptive neurons ([231]25) and an
   overall higher percentage of CGRP-positive peptidergic axons ([232]25).
   In addition, increased local IL-6 can, directly and indirectly,
   modulate responses in surrounding tissues, sensitizing the segmental
   DRGs and further perpetuating pain ([233]68). In both acute and chronic
   pain models, ECM modulation in the DRGs correlates with pain ([234]70).
   MMPs, including MMP2 and MMP9, which break down ECM, can interact with
   cytokines and integrins, reducing the local capacity for peripheral
   nerve regeneration, neuroplasticity, and connectivity and potentially
   contributing to neuropathy and chronic pain ([235]71). Moreover,
   astrocytes interact with the ECM by expressing several matrix
   components, such as hyaluronan, or proliferate in response to local
   fibronectin signaling ([236]72, [237]73). In agreement with this idea,
   we observed a concomitant modulation of astrocyte and matrix pathways
   in SM/J DRGs, with up-regulated DEGs linked to microglial migration,
   perception of pain signaling, regulation of neuron death, and
   neuroinflammatory response pathways.

   The spinal cord conducts pain-stimulated signals through the central
   nervous system to the regions of the brain ([238]74). We observed a
   higher abundance of IB4+ axons, a nociceptive neuronal marker, in the
   dorsal horn of SM/J. In addition, we noted increased GFAP and changes
   in astrocyte distribution in the deeper layers of the tissue. While
   higher levels of GFAP are associated with increased pain, it was
   unexpected to observe an age-associated change in the number and
   distribution of astrocytes in the dorsal horn in both BL/6 and SM/J
   strains ([239]31). Expression of pain-associated markers was restricted
   to the lumbar region of the spinal cord, with no differences observed
   in the thoracic region. Coalescence of disc herniation, DRG, and spinal
   cord changes in the lumbar spine of aged SM/J mice supports the
   hypothesis of a discogenic pain source in these animals.

   Transcriptomic studies of the DRGs additionally highlighted the
   associations between discogenic sensitization, inflammatory cytokine
   modulation, and local immune system regulation, associations that are
   well accepted ([240]75). We explored the cross-talk between changes in
   the DRGs and systemic inflammation, as aging and senescence are
   associated with increased systemic inflammation, reduced immune
   activity and efficiency, and increased prevalence of musculoskeletal
   diseases and neuropathic pain ([241]67, [242]76, [243]77). Several
   cytokines, such as IL-1, IL-10, IFN-γ, and TNF-α, have been associated
   with disc pathology ([244]78), with TNF-α and IL-1 playing prominent
   roles in disc herniation and homeostasis ([245]79). Our results showed
   a negative correlation between mechanical allodynia and systemic IFN-γ,
   IL-1, IL-2, and IL-4 levels. Considering that an insufficient and
   distorted immune response can limit hernia resorption and healing,
   these results align well with the histological evidence of unresolved
   herniations in aged SM/J mice. Noteworthily, local changes in several
   interleukins (IL-1, IL-2, IL-4, IL-6, IL-8, IL-10, IL-12, and IL-17) as
   well as IFN-γ regulation were observed in DRGs from aged SM/J mice,
   furthering the parallels between systemic and local inflammatory
   cytokines, immune regulation, and discogenic sensitization ([246]75).
   In addition, we found that in the AF and vertebral tissues, T cell
   activation correlated with persistent hernia. A previous study of
   herniated human discs described the presence of activated T cells in
   17% of the herniated discs ([247]80). Building on this framework, we
   analyzed immune cell profiles in the primary lymphoid organ (spleen)
   and peripheral circulation, which showed a contraction of the lymphoid
   compartment at the expense of the myeloid compartment during BL/6
   aging, as reported previously ([248]81). Our analyses showed an
   overactivation of CD8+ T cells and inhibition of CD4+ T cells in SM/J
   mice. These results agree with recent studies attributing a key role to
   CD4+ subsets, namely regulatory T cell and T[H]2, in suppressing pain
   and CD8+ T cells in chronic pain ([249]82, [250]83)[.] PBMC analysis
   suggests a genetic background–based dysregulation of CD4+ and CD8+
   cells, as SM/J mice show a unique population of these cell types by 2
   years of age. Unexpectedly, while we observed a higher expression of
   INF-γ in splenic CD8+ T cells, systemic levels showed an inverse
   correlation with pain. While active CD8+ T cells and myeloid cells
   express higher INF-γ in the spleen, the systemic levels likely reflect
   the total decrease in INF-γ in splenic Lin− cells and PBMCs, which
   showed lower levels of monocytes, NK cells, and T cells. Despite a
   lower proportion of T cells systemically, GSEA suggested an elevated T
   cell–driven response in 20-M SM/J, which aligns with findings in
   patients with chronic pain ([251]83, [252]84). In line with this, SM/J
   splenocytes showed as immune response skewing toward T[H]1, supported
   by lower levels of IL-2 and IL-4 ([253]82–[254]84). Sciatic nerve
   injury and osteoarthritis pain phenotypes are modulated by B cell
   recruitment and B cells expressing Mu opioid receptors ([255]85,
   [256]86); notably, SM/J mice also evidence higher proportions of B
   cells in circulation during aging, and these cells may contribute to
   pain sensitization associated with herniations. As previous studies
   reported a vital role of innate immune system quiescence in
   perpetuating herniation-associated back pain, we observed a significant
   deficiency of innate response, specifically the monocytes,
   granulocytes, and NK cells in the peripheral circulation of SM/J
   ([257]13). While we observed an increase of peripheral macrophages in
   the 20-M SM/J group with herniations, transcriptomic profiles showed a
   suppressed immune response, suggesting an insufficient immune response
   and possible sterile inflammation. A recent study specifically
   investigated disc hernia regression in the context of macrophage
   depletion and restoration, providing critical evidence of the need for
   sufficient macrophage activity in resorption ([258]12). The possibility
   of macrophage deficiency promoting sterile inflammation was further
   supported by our plasma analysis, showing decreased IFN-γ and lower
   intracellular cytokine abundance in the splenic innate immune cell
   populations analyzed with CyTOF. Acknowledging the differences in
   immune regulation between humans and mice, Moss et al. ([259]87)
   described a common aging-related macrophage defect that impairs
   migration and phagocytosis in both species. This finding further
   supports the concept of immune quiescence and may help explain the
   higher incidence of disc hernia in the aging population across species.
   Noteworthily, local increases of IL-6, in addition to systemic
   decreases of IL-2 and IL-4, are indicative of M1 macrophage
   polarization, which has been linked to chronic pain ([260]88). A
   previous study has found varying age-related trends in NK cell
   abundance, hypothesizing that this may reflect the lack of universality
   in NK cell dynamics during aging ([261]89). In the present study, CyTOF
   and PBMC analyses showed an age-associated reduction of NK cells in
   SM/J mice, underscoring an additional indicator of dampened innate
   immunity. Together, these findings highlight a possible cross-talk
   between increased CD8+ T cell activity and dampening of the innate
   immune system, driven by a dysregulated T[H]1 immune response, failed
   hernia resorption, and consequently, the resolution of back pain.
   Future preclinical and clinical studies will be critical in evaluating
   the relationship between immune signatures and disc herniation.

   In summary, our studies describe a novel clinically relevant mouse
   model of spontaneous, age-associated disc herniation on a wild-type
   strain background. As in humans, genetic susceptibility plays a
   profound role in SM/J disc pathology, and SM/J mice evidence several
   salient features of discogenic chronic pain ([262]Fig. 8A). Despite the
   observational nature, we can consider SM/J genetic background and age
   as independent variables and upstream drivers of the multisystem
   changes observed in the disc, nervous tissue, and immune system. Our
   studies provide new insights into the molecular pathophysiology of
   chronic discogenic pain and demonstrate a key role of immune cell
   dysregulation in both lymphoid organs and peripheral circulation in
   mediating, and as a predisposing factor to, this process.

METHODS

Mice, treatment, and study design

   All animal experiments were performed under Institutional Animal Care
   and Use Committee protocol 01679 approved by Thomas Jefferson
   University. Mice were studied at 6 (healthy adult, 6 M), 12
   (middle-aged, 12 M), and 20 M of age (aged, 20 M; mice were collected
   between 18 and 22 M, with an average age of 20 M). C57BL/6 mice were
   obtained from the rodent colony at the National Institutes of Aging.
   The SM/J mice were obtained from the Jackson Laboratory (stock no.
   000687, Jackson Labs), bred, and aged at Thomas Jefferson University.
   The terminal point of 18 to 22 M was chosen on the basis of the
   comparative life phases of C57BL/6J mice and humans ([263]90).

Histological analysis

   Histological analysis was performed using six lumbar levels from L1-S1:
   n[BL/6 6 M] = 5 (3 females and 2 males); n[BL/6 12 M] = 6 (2 females
   and 4 males); n[BL/6 20 M] = 9 (3 females and 6 males); n[SM/J 6 M] = 7
   (4 females and 3 males); n[SM/J 12 M] = 10 (6 females and 4 males);
   n[SM/J 20 M] = 8 (2 females and 6 males). At terminal time points, mice
   were either euthanized with CO[2] asphyxiation or anesthetized with
   ketamine followed by transcardiac perfusion with 0.9% saline, followed
   by 4% paraformaldehyde (PFA). Lumbar spines were dissected and fixed in
   4% PFA in phosphate-buffered saline (PBS) for 6 to 48 hours,
   decalcified in 20% EDTA, and embedded in optimal cutting temperature
   compound (OCT) or paraffin. Histological assessment was performed using
   10-μm frozen or 7-μm paraffin-embedded midcoronal sections from six
   lumbar levels (L1-S1). The spinal cord (thoracic and lumbar) and DRGs
   (lumbar) were harvested from PFA-perfused animals, cryoprotected in 30%
   sucrose for 3 days, and embedded in OCT. These tissue blocks were cut
   serially in the sagittal planes at a thickness of 30 μm. Sections were
   collected on glass slides and stored at −20°C until analysis.
   Safranin-O/fast green/hematoxylin staining was conducted on six lumbar
   levels (L1-S1) of each spine. Picrosirius red staining was used to
   visualize the collagen fibers within the disc. Staining was visualized
   using an Axio Imager 2 microscope (Carl Zeiss) using 5×/0.15
   N-Achroplan or 20×/0.5 EC Plan-Neofluar objectives (Carl Zeiss) or a
   polarizing light microscope (Eclipse LV100 POL; Nikon, Tokyo, Japan).
   Four blinded observers scored midcoronal sections from six lumbar discs
   per mouse using a modified Thompson grading scale ([264]19).

Immunohistology and cell number measurements

   Following antigen retrieval, deparaffinized sections were blocked in 5%
   normal serum in PBS with Tween 20 and incubated with antibodies against
   COL2 (1:400, Fitzgerald, 70R-CR008), COMP (1:200, Abcam, ab231977), CS
   (1:300, Abcam, ab11570), CA3 (1:150, Santa Cruz), IL-6 (1:50, Novus,
   NB600-1131), and MMP13 (1:150, Abcam, ab39012). A MOM kit (Vector
   Laboratories, BMK-2202) was used for blocking and primary antibody
   incubation for ARGxx (1:200, Abcam, ab3773) staining. Tissue sections
   were washed and incubated with species-appropriate Alexa Fluor
   594–conjugated secondary antibodies (Jackson ImmunoResearch Lab Inc.,
   1:700). Sections were mounted with ProLong Diamond Antifade Mountant
   with 4′,6-diamidino-2-phenylindole (DAPI; Thermo Fisher Scientific,
   [265]P36971) and visualized with an Axio Imager 2 microscope using
   5×/0.15 N-Achroplan or 20×/0.5 EC Plan-Neofluar objectives, and images
   were captured with an Axiocam MRm monochromatic camera (Carl Zeiss).
   Staining area and cell number quantification was performed using ImageJ
   software version 1.53e (last access: 6 September 2020;
   [266]http://rsb.info.nih.gov/ij/) ([267]91).

   Slides with spinal cord sections were first warmed for 1 hour at 37°C
   before three 5-min washes in PBS to remove OCT medium. Slides were then
   blocked for 1 hour in blocking buffer (5% Donkey Serum and 0.1% Triton
   X-100 in PBS) and then incubated at 4°C overnight with primary
   antibodies against Iba1 (1:600, WAKO, AB_839504) and GFAP (1:400, Dako,
   Carpinteria, CA; AB_10013482) diluted in blocking buffer. Slides were
   then washed and incubated in the appropriate secondary antibody (1:200,
   Donkey anti-Rabbit Alexa Fluor 647, AB_2492288, or Donkey anti-rabbit
   Alexa Fluor 488, AB_2313584, in blocking buffer) before being mounted
   with Fluorsave (Millipore Sigma) and imaged.

Neuronal, vascular, and immune staining in the spinal column

   Endomucin (1:300, Santa Cruz Biotechnology, sc-65495), neurofilament
   (1:200, BioLegend, 837904), IB4 (1:200, Sigma-Aldrich, L2895), CGRP
   (1:50, Abcam, ab36001), IB4 (1:200, Sigma-Aldrich, L2895), CD19 (1:50,
   BioLegend, 115569), CD3 (1:50, BioLegend, 100202), Ly6G (1:50,
   BioLegend, 127602), MHCII (1:50, Thermo Fisher Scientific,
   50-112-9377), and F4/80 (1:50, Cell Signaling, 70076S) stains were
   imaged on a Zeiss LSM 800 equipped with a Plan Apochromat 20×,
   0.8–numerical aperture air objective (Zeiss). Wavelength excitation was
   488, 594, 647, or 405 nm based on antibody conjugates, and samples were
   captured as tiled or untiled z-stacks at a 1.0-μm interval. Tiled
   images were oriented to an upright angle and cropped to 1024 by 1024,
   and brightness and contrast were linearly adjusted using ImageJ 1.52i
   (National Institute of Health).

TUNEL assay

   Cell viability staining was performed using the in situ cell death
   detection kit (Roche Diagnostic). In summary, sections were
   deparaffinized and permeabilized using proteinase K (20 μg/ml) for 15
   min, and the TUNEL assay was carried out per the manufacturer’s
   protocol. Sections were washed and mounted with ProLong Diamond
   Antifade Mountant with DAPI and visualized and imaged with the Axio
   Imager 2 microscope.

Digital image analysis

   All immunohistochemical quantification was conducted in grayscale using
   the Fiji package of ImageJ ([268]91). Images underwent thresholding to
   create binary images, and NP, AF, end-plate (EP), or vertebral
   compartments were manually defined using the Freehand Tool. These
   defined regions of interest were then analyzed using either the Analyze
   Particles (TUNEL and cell number quantification) function or the Area
   Fraction measurement.

Tissue RNA isolation and microarray analysis

   AF tissues were dissected from lumbar discs (L1-S1) of 12- and 20-M-old
   BL/6 and SM/J animals [n[BL/6 12 M] = 5 (5 males); n[BL/6 20 M] = 4 (3
   females and 1 male); n[SM/J 12 M] = 4 (4 males); n[SM/J 20 M] = 4 (2
   females and 2 males)]. Eight to 10 DRGs were collected from the lumbar
   region (L1-S1) from 12- and 20-M BL/6 and SM/J mice (three males and
   three females per strain per time point). Pooled tissue from a single
   animal served as an individual sample. Samples were homogenized, and
   total DNA-free RNA was extracted using the RNeasy Mini kit (Qiagen).

   Total RNA with an RNA integrity number >4 was used for microarray
   analysis. Fragmented biotin-labeled cDNA was synthesized using the
   GeneChip WT Plus kit (Thermo Fisher Scientific). Gene chips (Mouse
   Clariom S) were hybridized with biotin-labeled cDNA. Arrays were washed
   and stained with the GeneChip hybridization wash & stain kit and
   scanned on an Affymetrix Gene Chip Scanner 3000 7G. Quality control of
   the experiment was performed in Expression Console Software version
   1.4.1.. CHP files were generated by sst-rma normalization from
   Affymetrix.CEL files using Expression Console Software. Only
   protein-coding genes were included in the analyses. Detection above
   background higher than 50% was used for significance analysis of
   microarrays, and the FDR value was set at 5% and FC ≥ 2. Biological
   process enrichment analysis was performed using the PANTHER
   Overrepresentation Test with Panther GO Slim Biological Processes and
   Panther Pathway database annotations and a binomial statistical test
   with FDR ≤ 0.05. Analyses and visualizations were conducted in
   Affymetrix Transcriptome Analysis Console 4.0 software. Pathway
   schematic analyses were performed using the Transcriptome Analysis
   Console. Array data are deposited in the Gene Expression Omnibus (GEO)
   database ([269]GSE252383).

Human microarray analysis

   Human AF transcriptomic analysis was performed with the [270]GSE70362
   dataset available in the GEO database by Kazezian et al. ([271]23).
   First, degenerative and healthy clusters were obtained by hierarchical
   clustering of DEGs using Euclidean metrics. Human AF samples with
   Thompson grades I and II were grouped into a healthy category, those
   with grades III and IV were considered degenerated, and those with
   grade V were considered herniated. To maintain consistency between
   human and mouse samples, healthy versus degenerated was compared with
   1-year SM/J and BL/6 samples and herniated versus healthy and
   degenerated compared against SM/J 20-M aging DEGs, obtained with
   P ≤ 0.05 and FC ≥ 1.3 for healthy versus degeneration comparison and
   P ≤ 0.05 for herniation versus healthy and degeneration. For human and
   mouse comparison, both human and mouse DEGs were first mapped using
   Panther generic mapping from reference proteome genome organisms and
   then overlapped to identify the common DEGs between each group. To
   study the biological pathways commonly modulated between humans and
   mice (independent of specific genes), we used theme/concept analysis
   using the CompBio tool (PercayAI Inc., St. Louis, MO) (P ≤ 0.05; table
   S2) ([272]92).

Plasma cytokine analysis

   Blood was collected by intracardiac puncture in heparin-coated syringes
   and centrifuged at 1500 relative centrifugal force (rcf) at 4°C for 15
   min to isolate the plasma. Samples were stored at −80°C until the time
   of analysis. According to the manufacturer’s specifications, levels of
   pro-inflammatory proteins, cytokines, and T[H]17 mediators were
   analyzed using the V-PLEX Mouse Cytokine 29-Plex Kit (Meso Scale
   Diagnostics, Rockville, MD).

Behavior testing

   Modified versions of the Hargreaves, cold stimulus, and von Frey
   filament tests were used to measure heat, cold, and mechanical
   sensitivity, respectively ([273]28, [274]29). Hind paw withdrawal for
   Hargreaves, cold, and von Frey tests was defined by removing the paw
   from the infrared, dry ice (through a glass plate), or mechanical
   stimulus. All scored withdrawal responses were associated with one or
   more supraspinal responses, including vocalization and licking and/or
   guarding of the tested paw. Spontaneous movements were not recorded,
   and each trial was conducted in triplicate. Cutoffs of 30 s for the
   Hargreaves and cold stimuli tests and 2 g of force for the von Frey
   test were used to prevent damage to the animal’s paw.

Subchondral bone microcomputed tomography analysis

   Lumbar spines were dissected from 20-M C57BL/6J and SM/J mice
   [n[BL/6] = 5 (2 females and 3 males); n[SM/J] = 5 (3 females and 2
   males)], fixed in 4% PFA for 48 hours, and rehydrated with 1× PBS for
   assessment with the SkyScan 1275 μCT imaging system (Bruker). Scanning
   settings were 200-μA current, 50-kV voltage, and 15 μm per pixel. The
   skeletal parameters assessed by microcomputed tomography followed
   previously published nomenclature guidelines ([275]93). Images were
   reconstructed through NRecon (Bruker) before analysis in CTAn (Bruker).
   The regions of interest were limited to subchondral bone and end plate
   by selecting five cross-sectional slices, which excluded extraneous
   trabeculae of the vertebral body. The density of mineralized tissue was
   obtained through CTAn using appropriate thresholding. Subsequent
   two-dimensional analysis was used to assess bone volume (BV), mean bone
   area of selected cross sections (B.Ar), cross-sectional thickness
   (Cs.Th), bone surface fraction (BS/BV), and closed porosity ([276]93).

PBMC single-cell preparation and sequencing

   Whole blood was collected from 12- and 20-M SM/J and BL/6 mice (n[12 M
   BL/6] = 2 females and 2 males; n[12 M SM/J] = 2 females and 1 male;
   n[20 M BL/6] = 2 females and 2 males; n[20 M SM/J] = 4 males) via
   intracardiac puncture and diluted into an equal volume of Dulbecco’s
   PBS and 2% fetal bovine serum (FBS; StemCell Technologies, 07905).
   According to the manufacturer’s specifications, blood was layered on
   top of Lymphoprep (StemCell Technologies, 07851) in a 15-ml tube and
   centrifuged at 800 rcf for 20 min at room temperature. The mononuclear
   cells were washed twice with PBS + 2% FBS, and platelets were removed
   with an additional centrifugation at 120 rcf for 10 min. To facilitate
   the removal of any remaining red blood cells (RBCs), a lysis step was
   performed (BioLegend, 420302). Briefly, lysis buffer (1×) was added to
   the cell suspension and incubated for 5 min, with periodic inversion.
   Cells were then centrifuged at 350 rcf for 10 min, and PBMCs were
   washed twice with PBS + 2% FBS. Cell suspensions were processed for
   scRNA-seq using the 10X Genomics 3’ version 3 kit (10x Genomics,
   Pleasanton, CA), with interstrain multiplexing, according to the
   manufacturer’s specifications. Sequencing was performed with Illumina
   NextSeq 2000 according to the manufacturer’s instructions (Illumina).

scRNA-seq data analysis

   Raw data were aligned to the mouse reference transcriptome
   (mm10-2020-A) and underwent preprocessing with Cell Ranger Multi
   version 7.0.1. According to the authors’ recommended workflow, the
   count matrix was analyzed using Seurat version 4.3.0 ([277]94).
   Briefly, data were converted into Seurat Objects using the
   CreateSeuratObject function, and all cells were analyzed for their
   unique molecular identifier and mitochondrial gene fractions; cells
   were included if 200 < n[feat] < 2500, and mitochondrial genes
   accounted for <5%. Data were integrated with logarithmic normalization
   (LogNormalize) and organized according to variation in the data using
   FindVariableFeatures. Dimensional reduction was performed using
   principal components analysis, and the FindNeighbors (dims = 1:15)
   function, followed by the FindClusters functions, was used to calculate
   the number of cell clusters using the Louvain algorithm, with the
   resolution parameter set to 0.5. DimPlot was then used to visualize the
   data as UMAP. Cell types for each age and strain condition were
   determined using SingleR, with the MouseRNAseqData reference dataset
   and main labels ([278]33). Individual cell types were individually
   subset and analyzed via the same methodology to determine strain- and
   age-based variation within each cell type, with DEG results in table
   S3. The clusterProfiler package was used to conduct GSEA on DEGs
   (FDR < 0.05) identified in various strain- and age-based comparisons
   ([279]36). Data are deposited and may be accessed in the GEO database
   (GSE252474), and relevant R scripts can be found at
   [280]https://github.com/mrisbudlab or
   [281]https://datadryad.org/stash/share/UIzO6e_x3uLqiw3zPr0Zv0ntHbXZPrXP
   -hG02aA5Iec.

Splenocyte preparation and CyToF

   Spleens were collected from 12- and 20-M SM/J and BL/6 mice (n[12 M
   BL/6] = 2 females and 2 males; n[12 M SM/J] = 2 females and 1 male;
   n[20 M BL/6] = 2 females and 2 males; n[20 M SM/J] = 4 males).
   Splenocytes were dissociated from the tissue matrix by gentle
   compression through a 100-μm strainer (StemCell Technologies, 27217).
   RBC lysis (BioLegend, 420302) was performed, followed by two PBS + 3%
   FBS washes. Cells were then resuspended in complete medium
   (Sigma-Aldrich, D6046) + 10% FBS (Sigma-Aldrich, F2442) + 2% 50×
   penicillin/streptomycin (Corning, 30-001-CI) at a concentration of 1 ×
   10^6 cells/ml. To stimulate intracellular cytokine production, Cell
   Activation Cocktail (2 μl/ml) with Brefeldin A (BioLegend, 423303) was
   incubated with the cells at 37°C for 4 hours. After incubation, the
   mixture was washed twice in serum-free media and resuspended in
   serum-free media at a concentration of 2 × 10^7 cells/ml. To label dead
   cells, Cell-ID Cisplatin (Standard Biotools, 201064) was added to each
   cell suspension at a concentration of 10 μM. Cell staining according to
   the manufacturer’s protocol for the Maxpar Mouse Sp/LN Phenotyping
   Panel Kit (16 markers: Ly6G/C[Gr1], CD11c, CD69, CD45, CD11B[MAC1],
   CD19, CD25, CD3e, TER-110, CD62L, CD8a, TCRI^2, NK1.1, CD44, CD4, and
   B220) (Standard Biotools, 201306) was carried out for all samples. 18-M
   SM/J and BL/6 splenocytes underwent a subsequent round of staining
   according to the manufacturer’s specs for the Maxpar Mouse
   Intracellular Cytokine I Panel Kit (eight markers: IFN-γ, IL-2, IL-4,
   IL-5, IL-6, IL-10, IL-17A, and TNF-α) (Standard Biotools, 201310).
   After all staining, samples were washed, fixed with 1.6% PFA in PBS,
   and incubated with Cell-ID Intercalator-Ir, per the Maxpar staining
   protocols. After this incubation, cells were washed once, passed
   through a 100-μm cell strainer, and resuspended in freezing media
   composed of 90% FBS + 10% dimethyl sulfoxide. Cells were stored at
   −80°C until the time of analysis. Cells were quickly thawed at 37°C and
   washed three times in CAS buffer at 800g for 10 min. Cells were
   resuspended at 1 × 10^6/ml in CAS buffer and strained through a 35-μm
   filter before acquisition. Each sample was acquired on a CyTOF 2 Helios
   mass cytometer.

CyToF analysis

   Single-cell data have been clustered using the FlowSOM R package and
   labeled using the Ek’Balam algorithm ([282]95, [283]96). Cell subset
   definitions follow guidelines by Maecker et al. ([284]97) and Finak
   et al. ([285]98). Cluster labeling, method implementation, and
   visualization were conducted through the Astrolabe Cytometry Platform
   (Astrolabe Diagnostics Inc.). The Astrolabe Mass Cytometry Platform is
   a cloud-based service for analyzing mass cytometry data. Astrolabe
   labels events into immune subsets in a two-stage process. First, events
   are clustered using self-organizing maps implemented by the FlowSOM
   package. Clusters are then labeled, following a manually curated gating
   hierarchy. Multidimensional scaling maps were generated using the
   cmdscale R function ([286]99). Differential abundance analysis was done
   using the edgeR package ([287]100) following the methods outlined by
   Lun et al. ([288]101). Differential expression analysis was done using
   the limma R package ([289]102) followed by diffcyt ([290]103).

Statistical analyses

   All statistical analyses were performed using Prism7 (GraphPad, La
   Jolla). Data distribution was assessed with the Shapiro-Wilk normality
   test, and the differences between the two groups were analyzed by t
   test or Mann-Whitney test, as appropriate. Differences between more
   than two groups were analyzed by analysis of variance (ANOVA) or
   Kruskal-Wallis test for non-normally distributed data, followed by
   Dunn’s multiple comparison test. Chi-square test was used to analyze
   the differences between the distributions of percentages. P ≤ 0.05 was
   considered a statistically significant difference.

Acknowledgments