Progressive abnormalities found in aged zebrafish vertebral columns. a 3D rendering from μCT images of young (1 year) and aged (3 years) spines. AV anterior vertebra, PV posterior vertebra, IVD (arrow). Scale bar = 500 μm. b Frequent changes observed at the IVD. 1: normal IVD; 2: osteophytes (pink arrow), vertebral misalignment and IVD narrowing (green arrow); 3: endplate sclerosis; 4: sclerosis (pink arrow) and IVD narrowing (green arrow); 5: sclerosis and fusion (pink arrow); 6: sclerosis and IVD calcification (pink arrow). Scale bar = 100 μm. c Heat map graph showing spinal morphological changes classified by severity during aging (1 year n = 36 (42% females, 58% males), 2 years n = 16 (57% females, 43% males), and 3 years n = 34 (55% females, 45% males)). Average fish standard lengths (measured from tip of the head to the last vertebral column): 1 year = 3.22 cm (0.19 SD), 2 years = 3.41 cm (0.18 SD), and 3 years = 3.46 cm (0.23 SD). d 3D rendering from μCT images of young and aged fish, color coded to show bone mineral density changes. The selected area of the spine (dashed box) is magnified, as shown on the right of the panel. Higher density colocalizes with regions of sclerosis and deformities (arrowhead) in the aged spine. Scale bar = 500 μm. e TMD (tissue mineral density) retrieved from the third thoracic vertebrae in young (1 year) and aged (3 years) fish. Nonparametric, two-tailed, Mann–Whitney test; data are the mean and SD. P values are indicated. f 3D volume rendering from a μCT image of wt fish showing an individual vertebral centrum segmented by computational automation. g Relative bone density from the vertebral centra in aging fish (3 months to 3 years). Average standard lengths: 3 months = 2.6 cm (0.18 SD); 3 years = 3.46 cm (0.23 SD). The notch plot was scaled by the average value from the 3-year centra. h Within-sample standard deviation in bone density. The notch plot was scaled by its average value from the 3-year centra. c, e Generated in Prism 8. g, h Graph was generated in Python

Altered vertebral column biomechanics in aged and sp7^−/− zebrafish show increased bone fragility. a Intact vertebral column motion segments (composed of three vertebrae surrounded by trunk musculature) were placed in a material testing stage (MTS), followed by μCT imaging at increasing compression forces. An example of vertebral column segments before compression (N = 0) and after failure (N > 0) is shown. Anterior and posterior orientations of the zebrafish are annotated. b 3D volumetric rendering from μCT images of 3-month-old wt and sp7^−/− zebrafish color coded for tissue mineral density (TMD). The dashed box shows the magnified region presented on the right of the panel. Scale bars = 500 μm. c Volumetric rendering from μCT images of 1-year-old wt and sp7^−/− zebrafish. Note IVD calcification (arrowhead) and abnormal shapes of the centra (arrow) in sp7^−/− zebrafish. Scale bars = 500 μm. d Failure points during compressive forces of 3-month-old (3 mo) sp7^−/− zebrafish, wt zebrafish siblings, 1-year-old wt siblings, and 3-year-old wt siblings (n = 3 per group). One-way ANOVA and post hoc Holm–Sidak’s multiple comparisons test were used; data are the mean and SD. P values are indicated. The graph was generated in Prism 8. Scale bars = 500 μm. e TMD of 3-month-old (3 mo) wt and sp7^−/− zebrafish (wt n = 5 (3 males +2 females), sp7^−/− zebrafish n = 3 (2 males + 1 female)). Standard fish lengths: wt = 2.43 cm (0.2 SD); sp7^−/− = 2.17 (0.4 SD). Nonparametric, two-tailed, Mann–Whitney test; data are the mean and SD. P values are indicated. f Radiograph of a 16-year-old male with a frameshift mutation in SP7 (c.1052delA). Frontal and lateral images are displayed. Regions 1 and 2 are magnified and displayed on the right. Lumbar spinal curvature (asterisk), biconcave vertebrae (pink arrows), the spread of IVD, and signs of IVD calcification (cyan arrowhead) are shown

SRCT reveals subtle bone morphological abnormalities and alterations in mineral density distribution. a 3D volumetric rendering from μCT images showing fine endplate abnormalities. 1: normal IVD; 2: endplate sclerosis (blue arrow) and IVD calcification (pink arrow); 3: sclerosis (blue arrow), uneven edges and AF calcification (pink arrow); 4: sclerosis (blue arrow), AF calcification (pink arrow); osteophyte (dashed blue arrow). Scale bars = 50 μm. b 3D volumetric rendering from μCT images showing IVD internal changes in aged and sp7^−/− zebrafish. Calcification depicted by pink arrow (SRCT 3D rendering). Note, young zebrafish with IVD and disc misalignment shows points of IVD calcification. Virtual sections (orthoslices) color coded for bone density (SRCT orthoslice). Dashed box regions are magnified and displayed on the right. Scale bars = 50 μm. c Quantification of mineralization robustness (%) in wt young (1 y), wt aged (3 y), and sp7^−/− 1 y zebrafish (n = 3 per group, average from 4 sections per fish). One-way ANOVA and post hoc Holm–Sidak’s multiple comparisons test were used; data are the mean and SD. P values are indicated. d Paraffin sections stained with picrosirius red showing the thickness of collagen fibers (thin = blue/green and thick= red/orange). Scale bars = 50 μm. e Quantification of collagen fiber thickness from picrosirius red staining (blue = thin, green = medium, red = thick) (1 y, n = 8; 3 y, n = 8; 3 months, n = 7; sp7^−/− 1 y, n = 6; we analyzed three sections per fish). One-way ANOVA, post hoc Tukey’s multiple comparisons test; data are the mean and SD. P values are indicated

The osteocyte lacunar profile is unchanged in aged fishbut is dramatically compromised in young sp7^−/− fish. a The top panel presents 3D volumetric rendering and postimaging analysis of young (1 y wt), aged (3 y wt), and sp7^−/− (1 y) (lateral view of the IVD, 0°), showing vertebral bone and segmented lacunae (labeled in green), resulting from automated image segmentation using deep learning. Note the distribution of lacunae at the endplates and along the bone, with a dramatic reduction in sp7^−/−. The bottom panel shows segmented lacunae (90° clockwise rotation from the lateral view of the IVD), color coded to show orientation (angle), with the center of the vertebral centrum as reference for orientation. Scale bars = 50 μm. b Violin plot of lacunae orientation from the center of the centrum. Data are individual zebrafish (n = 3). Nested one-way ANOVA, post hoc Tukey’s multiple comparisons test. P values are indicated. ns nonsignificant. c Violin plot of number of lacunae per volume of bone. One-way ANOVA, post hoc Tukey’s multiple comparisons test. P values are indicated. ns nonsignificant. d Violin plot of lacunae distances from the bone surface. Nested one-way ANOVA, post hoc Tukey’s multiple comparisons test. P values are indicated. ns nonsignificant. e Violin plot of the distance to the closest lacunae. Data are the mean (n = 3). Nested one-way ANOVA, post hoc Tukey’s multiple comparisons test. P values are indicated. ns nonsignificant. f Violin plot of lacunae volume. Note the subpopulation of lacunae showing a small volume in sp7^−/− zebrafish. Data are the mean (n = 3). Nested one-way ANOVA, post hoc Tukey’s multiple comparisons test. P values are indicated. ns nonsignificant. g The top 25% smallest lacunae volume from young (1 y), aged (3 y), and sp7^−/− (1 y) zebrafish were compared. One-way ANOVA, post hoc Tukey’s multiple comparisons test. P values are indicated. ns nonsignificant. h Violin plot of lacunae sphericity (circularity). Data are the mean (n = 3). Nested one-way ANOVA, post hoc Tukey’s multiple comparisons test. P values are indicated. ns nonsignificant

IVDD histopathology underlying 3D disc changes in zebrafish. a Schematic of the vertebral segments of the zebrafish. IVDD histopathology underlying 3D disc changes in zebrafish. The centrum and IVD regions are annotated. bo, bone; cf, collagen type I fiber bundles; dc, dense collagen type I matrix; el, elastin; ev, extracellular vacuole; fb, fibroblasts; nc, notochord cells, vacuolated cells and notochord epithelial cells attached to the notochord sheath; ne, notochord epithelium; ns, notochord sheath; ob, osteoblasts; se, septum; st, notochord strand. b Toluidine blue, pankeratin and pan cadherin immunostaining of young (1 y), aged (3 y), and sp7^−/− (1 y) discs. Toluidine blue of a middle section of the disc showing the NP in young, aged, and sp7^−/− discs (n = 3 per group, serial sections were analyzed). Note the fibrous and disorganized NP in aged and sp7^−/− (magenta arrows) discs. Higher magnification, with a focus on the AF, is shown on the right. Note collagen bundle fibers disrupted in the aged AF and loose organization in the sp7^−/− disc (magenta arrows). Pankeratin showing bone (yellow), keratin (orange), and glycosaminoglycans (green). Note the accumulation of keratin in aged discs and dehydration (weak stain for glycosaminoglycans) in aged and sp7^−/− discs. Immunostaining for pan cadherin showing the loss of cellularity in aged discs (arrow) and disorganized NPs in sp7^−/− discs, with vacuolated cells (pink arrowhead), acellular regions (pink arrow), and cellular agglomerates (pink dashed arrow) (n = 3 per group, ≥3 sections analyzed). Scale bars = 50 μm. c 3D volumetric rendering from enhanced (E) contrast μCT images of wt (1 y) and sp7^−/− (1 y) zebrafish. Note fibrous tissue organization in sp7^−/− zebrafish. Arrows point to an abnormal notochord sheath layer. Scale bars = 50 μm

Altered collagen fiber organization and disc herniation in zebrafish IVDD. a Max projection of image stacks from second harmonic generation (SHG) imaging color coded by depth. Individual collagen fiber was traced using CTFire. 3D volumetric rendering from SHG images were used to visualize disc herniation in 3D using Amira (gray = back-scatter SHG; purple = forward-scatter SHG). Orientation of the vertebral column is shown at the top left of the panel. The annuls fibrosus (AF) is indicated with a dashed line. Note the variation in AF length among the groups. Smooth perpendicular fiber organization at the endplate was observed only in young samples (white arrow, SHG). Thicker fibers were observed in the AF of aged and sp7^−/− zebrafish (dashed arrow). Bulging discs were found in aged zebrafish and severely bulging discs were found in sp7^−/− zebrafish (green arrows, 3D render). Scale bars = 50 μm. b Relative frequency distribution (%) of angles of collagen fibers from each group. The graph was generated in Prism 8. c Angle of the most common collagen fibers (75% of total fibers) show abnormal orientation in sp7^−/− fish (n = 3 fish per group). One-way ANOVA, post hoc Tukey’s multiple comparisons test; data are the mean and SD. P values are indicated when significant (P < 0.05). d Scanning electron microscopy (SEM) of 1-year-old wt and sp7^−/− zebrafish (1 y). Note a smooth transition from the vertebral bone to the AF in wt zebrafish and a cliff in sp7^−/− zebrafish (arrow). The orientation of the vertebral column is shown on the top left. Scale bars = 50 μm

Increased bone density in cathepsin K crp accompanied by premature IVDD. A Radiographs and 3D volumetric rendering from μCT images of 1-year-old wt zebrafish and 1-year-old ctsk crps, color coded for bone density (TMD). Note the dramatic IVD calcification in ctsk crps (arrow). Scale bars = 500 μm. B TMD showing differential density in ctsk crps [wt n = 6 (3 males, 3 females), ctsk crpsn = 8 (4 males, 4 females)]. Standard fish lengths: wt = 3.5 (0.3 SD); ctsk crp = 3.5 cm (0.09 SD). The graph was generated in Prism 8. Nonparametric, two-tailed, Mann–Whitney test; data are the mean and SD. P values are indicated. C Relative bone density from the vertebral centra in aging fish. The notch plot was scaled by the average value from the 3-year-old centra (cyan = 3 months wt; green = 1-year old wt; yellow = 1-year-old ctsk crps; orange = 2-year-old wt; magenta = 3-year-old wt). The graph was generated in Python. D Within-sample standard deviation if bone density. The notch plot was scaled by its average value from the 3-year-old centra (cyan = 3-month-old wt; green = 1-year-old wt; yellow = 1-year-old ctsk crp; orange = 2-year-old wt; magenta = 3-year-old wt). The graph was generated in Python. E Cross-correlation between mineral density and wt morphology of the discs was calculated for wt, ctsk crp and wt standard samples. Pearson correlation coefficient = −0.85, P = 7.8e⁻¹⁴. The graph was generated in Python. F Proposed U-shaped model to describe the association between BMD and IVDD. Either an abnormal increase or decrease in BMD is linked to IVDD. G Failure point during vertebral compression. Nonparametric, two-tailed, t-tests; data are the mean and SD. P values are indicated. The graph was generated in Prism 8. H Histological sections of wt and ctsk crp (1 year) zebrafish stained with toluidine blue and picrosirius red, with a focus on the NP and AF regions. Higher magnification of the AF from another histological section is shown on the right. Disorganized NP (dashed arrows, toluidine blue) and abnormal AF, displaying loss of the elastin layer (arrows, toluidine blue), were detected in the affected IVD of ctsk crps. Picrosirius red staining shows different colors in the bone of ctsk crps (higher magnification), indicating bone quality impairment (white arrow). Scale bars = 50 μm. I Quantification of the thickness of collagen fibers from picrosirius red staining (blue = thin, green = medium, red = thick) (1-year-old n = 8, ctsk crpn = 8). Parametric, two-tailed, T-test; data are the mean and SD. P values are indicated. The graph was generated in Prism 8