Muscular loading affects the 3D structure of both the mineralized rudiment and growth plate at early stages of bone formation

Fetal immobilization affects skeletal development and can lead to severe malformations. Still, how mechanical load affects embryonic bone formation is not fully elucidated. This study combines mechanobiology, image analysis and developmental biology, to investigate the structural effects of muscular loading on embryonic long bones. We present a novel approach involving a semi-automatic workflow, to study the spatial and temporal evolutions of both hard and soft tissues in 3D without any contrast agent at micrometrical resolution. Using high- resolution phase-contrast-enhanced X-ray synchrotron microtomography, we compare the humeri of Splotch-delayed embryonic mice lacking skeletal muscles with healthy littermates. The effects of skeletal muscles on bone formation was studied from the first stages of mineral deposition (Theiler Stages 23 and 24) to just before birth (Theiler Stage 27). The results show that muscle activity affects both growth plate and mineralized regions, especially during early embryonic development. When skeletal muscles were absent, there was reduced mineralization, altered tuberosity size and location, and, at early embryonic stages, decreased chondrocyte density, size and elongation compared to littermate controls. The proposed workflow enhances our understanding of mechanobiology of early bone formation and could be implemented for the study of other complex biological tissues.


Introduction
Fetal immobilization can have dramatic effects on the skeleton [1,2]. Reduced intrauterine movements can lead to smaller, thinner and weaker long bones that in some cases even fracture before birth, as well as to craniofacial and limb deformities and abnormal joint contractures [3][4][5]. Bone is believed to be deposited where it is mechanically required, i.e. when mechanical stimuli exceed a certain threshold [6,7]. Consequently, in the absence of muscular stimuli, the normal specific pattern of mineral deposition in bones can be altered [1]. Genetically modified mice with absent or non-contractile musculature are suitable to study the specific effects of an altered mechanical environment on mammalian bone development [8][9][10][11][12][13]. It has been shown that long bones of mice lacking functional muscles were significantly less mineralized than control bones and showed an irregular morphology [10]. This raises interesting questions such as: does mechanical contraction affect not only the bone morphology overall but also bone internal microstructure, and how is the primary bone organized in the bone volume when muscular contractions have been absent?
In mice, different long bones and joints seem to be affected differently by the lack of skeletal muscles [8]. This may be because the socalled "muscleless limb" mouse embryos still experience some passive mechanical stimuli, coming from movements of the mother or normal littermates [14]. Passive movements induce different levels of stimuli to the hind and forelimbs, potentially contributing to the humerus being more affected than the femur in muscleless limb mice [14]. However, at early embryonic stages the forces exerted by movements of the littermates are likely to be less impactful than later in development. An interesting question is whether mechanical forces are of equal importance throughout the stages of bone formation? Long bones grow through endochondral ossification, which is a process that is directed at the growth plates. The chondrocytes in the cartilage template progress from a resting state into a proliferating state, which is followed by hypertrophy and subsequent mineralization [4]. Sharir et al. describe how the embryonic bone growth starts with the deposition of a cylindrical stratum of mineral around the cartilage center in mice femora [2]. They show that the bone collar expands centrifugally, in a sequential manner, and that also bone resorption starts early. Their study reveals that osteoblasts respond to mechanical contractions and they propose that this could be the mechanism by which muscular loading regulates bone morphology [2]. This important study sets the foundation for further investigations, for instance, to understand if muscular contractions affect long bones uniformly, or if their effects are directional or localized (i.e. more pronounced along the bone longitudinal or transverse planes, or more mineralization in specific locations). Furthermore, a study from Stern et al. shows that long bones grow by isometric scaling [15]. During isometric bone growth, superstructures (such as deltoid tuberosity and condyles) maintain their relative positions, which require a constant balance between the proximal and distal growth rates [15]. It is however yet to be determined how lack of muscles affects the growth balance at the growth plates and if the relative positions of the deltoid tuberosity and condyles are maintained even when mechanical loading is absent. Therefore, it is important to study not only the already mineralized tissue but also how the forces affect the pre-mineralized tissue, i.e. the chondrocytes in the primary growth plates of the long bones. In recent years, some studies were conducted to understand the relationship between mechanical loading and endochondral ossification focusing on chondrocytes. Killion et al. concluded that in paralyzed young mice there are fewer and shorter chondrocyte columns [16]. Schwartz et al. hypothesized that muscle contractions are necessary in embryos to trigger normal chondrocyte formation as well as their alignment into columns, and that the lack of orientation could prevent column intercalation [17].
Despite these recent studies, the effect of mechanical loading on the chondrocyte properties in the primary growth plates are still not completely elucidated. Some studies have imaged cartilage in 3D after staining using optical projection tomography, confocal microscopy and micro X-ray computed tomography [18][19][20]. However, 3D imaging techniques suitable to address questions on both mineralized and cartilaginous tissue have only recently emerged, and many studies were, until now, based on 2D histology. Synchrotron X-ray tomography is a non-destructive method for 3D visualization and quantitative analysis of materials. In conventional X-ray tomography the contrast is given by the materials X-ray absorbance, thus imaging soft biological materials is very challenging. Phase-contrast imaging was developed to image materials that weakly absorb X-ray radiation and/or have similar density to each other [21]. The contrast is based on the difference between refractive indexes of materials, where the phase shifts that occur at the boundaries leads to 'edge enhancement' [21]. The technique has successfully been applied to soft biological tissues, e.g. lungs, heart, breast, cartilage, vessels, and nerves [22][23][24][25][26][27]. However, to our knowledge it has not been used to investigate the structural effect of mechanical forces on embryonic early bone deposition and to study growth plates.
This study aims to characterize and quantify the effect of absent skeletal muscles on the early mineralization of embryonic long bones by using high resolution phase-contrast synchrotron X-ray tomography. We present a new, semi-automatic framework to investigate both the spatial (3D) and temporal evolutions of the mineralization. Specifically, we investigate how the 3D microstructure of both mineralized bone and cellular components in the growth plate are affected by lack of mechanical loading. We focus on the humerus, as humeri are more affected than femurs by the lack of muscular stimulation [14], and compare the forming humerus of Splotch-delayed embryos of mice lacking skeletal muscles (mutants), with the humerus of littermates with normal skeletal muscles (controls). Embryonic stages were selected to study the effect of lack of muscular contractions on the early Theiler Stage (TS) when mineral deposition initiates (TS23 and TS24) and just before birth (TS27).

Animals
Pax3 Sp-d/Sp-d (Splotch-delayed) mutant mice were interbred by spontaneous mating and offsprings were subsequently genotyped as described previously, where ~25% of the embryos are mutants [28]. The absence of the transcription factor Pax3 prevents muscle progenitor cell migration into the limb buds, which in turn prevents the development of skeletal muscles [29]. Embryos were harvested and staged using Theiler morphological criteria [30]. This study included a total of 26 samples, comparing mutant embryos at stages TS23 (n = 4), TS24 (n = 4) and TS27 (n = 5) with littermates with normal skeletal muscles (controls) at TS23 (n = 4), TS24 (n = 5) and TS27 (n = 4). The experimental protocol was conducted in accordance with European legislation (Directive 2010/63/EU), project license number P39D18B9C. The left forelimbs were carefully dissected, fixed in 4% paraformaldehyde, and dehydrated in increasing ethanol gradient. Limbs were stored in the freezer in absolute ethanol until imaging. Samples were mounted in Eppendorf tubes on CryoCaps (Molcular Dimensions Ltd., UK) modified to remove the central pillar.

Synchrotron Phase contrast X-Ray Tomography
First, 24 samples were imaged with Synchrotron Phase contrast X-Ray Tomography at the Diamond-Manchester Imaging Branchline I13-2 of the Diamond Light Source (DLS) synchrotron (Oxfordshire, United Kingdom) [31,32]. A partially-coherent, near-parallel, polychromatic 'pink' beam (circa 8-30 keV) was generated by an undulator in an electron storage ring of 3.0 GeV voltage and 300 mA current. For data collection, the undulator gap was set to 6.0 mm. The beam was reflected from the platinum stripe of a grazing-incidence focusing mirror and high-pass filtered with 1.3 mm pyrolytic graphite, 3.2 mm aluminium and 70 μm steel, resulting in a beam of weighted-mean photon energy of ~29 keV. Samples were aligned for data collection under low-dose conditions (~10 min per sample) by temporarily setting the undulator gap to 12 mm. Slits were used to crop the beam just outside the field of view; this limited both sample exposure and the intensity of noise arising from scintillator defects. Samples were placed at room temperature on a HUBER 1002 manual goniometer (HUBER Diffraktionstechnik GmbH & Co. KG, Rimsting, Germany) (adjusted to make samples approximately vertical), mounted on perpendicular Newport MFA-PPD (Newport Corp., Irvine, CA, USA) linear stages, atop an Aerotech ABRT-260 (Aerotech Inc., Pittsburgh, PA, USA) rotation stage. Various propagation distances were tested, and ~190 mm was chosen to give a sufficient level of phase contrast. 8001 projections were acquired at equallyspaced angles over 180 • of continuous rotation ('fly scan'), with an extra projection (not used for reconstructions) collected at 180 • to check for possible sample deformation, bulk movements and beam damage relative to the first (0 • ) projection. Projections were collected by a pco. edge 5.5 Camera Link (PCO AG, Kelheim, Germany) detector (sCMOS sensor of 2560 × 2160 pixels) mounted on a visible light microscope of variable magnification. Magnification was controlled via rotation of a turret incorporating various scintillator-coupled objective lenses. A 2× objective (used for TS27 embryos) coupled to a 500 μm CdWO 4 scintillator that was mounted ahead of a 2× lens, provided 4× total magnification with a field of view (FOV) of 4.2 × 3.5 mm 2 and an effective pixel size of 1.625 μm. A 4× objective (used for TS23 and TS24 embryos and TS27 growth plates) coupled to a 500 μm CdWO 4 scintillator provided 8× total magnification, a FOV of 2.1 × 1.8 mm 2 and an effective pixel size of 0.8125 μm. Exposure time was 100 ms for the 2× objective lens, leading to scanning times of circa 15 min per sample. For the 4× objective, 150 ms exposure time was required, leading to a scanning time of circa 25 min per sample.
Additionally, two samples (one TS23 mutant and one TS27 mutant), were imaged at the X02DA TOMCAT beamline at the Swiss Light Source (SLS), Paul Scherrer Institute (Villigen, Switzerland) [33]. A High Numerical Aperture Microscope (custom made objective (optiquepeter.co m)), with 4× magnification, a FOV of 4.2 × 3.5 mm 2 and an effective pixel size of 1.63 μm, was coupled to a scintillator screen for X-ray to visible light conversion. A propagation distance of 39 mm was used to achieve phase contrast. The beam energy was set to 20 keV, and 1501 projections were acquired over 180 • of continuous rotation. Exposure time was 320 ms, leading to scanning times of circa 10 min per sample. Projections were collected by a pco.edge 5.5 Camera (detector sCMOS sensor of 2560 × 2160 pixels).
Projection images were flat-and dark-field corrected, and ring artifact suppression was performed [34]. Standard Paganin filtered back projection is used to reconstruct the volumes [34]. An example of the reconstructed data is shown in Video S1.

Image analysis
Initial processing of the tomographic images was performed using ImageJ [35] and Dragonfly 3.6, ORS [36]. To reduce file sizes and make subsequent processing more manageable, tomograms were downsampled from 32 to 8 bit and cropped. First the long bone axis of each sample was aligned to the tomogram's z-axis (ImageJ and Matlab scripts), and then the image stacks were divided and re-saved as two substacks comprising 1) the mineralized rudiment and 2) the growth plate. The two sub-stacks were analyzed separately using in-house Matlab scripts. The sample group sizes used for the two analyses are shown in Table 1.

Mineralized rudiment
The workflow to obtain quantitative information such as diameters, lengths and bone content of all samples is described in Fig. 1 and consisted of: A) Bone segmentation: while the mineralized embryonic rudiment is not yet truly bone, we use the terminology of bone for brevity. In order to remove noise from the images and improve the precision of the segmentation, a 2D median filter was applied sequentially in all 3 planes using a disk of variable size for each sample as a morphological structuring element adapted to the noise level in each image. Then, two binary volumes of the mineralized parts of the tomographic images were generated (Fig. 1A). First, the mineralized primary bone was selected by contrast threshold (rendered in green in Video S2). Due to higher noise level in few Table 1 Sample group sizes used for the analyses of the mineralized rudiment and the growth plate. The samples used for the growth plate analysis were a subset of those used for mineralized rudiment analysis.  samples (one control TS23, two controls TS24, and three mutants TS24), minor manual adjustments to the segmentation were performed using Dragonfly 3.6, ORS software. Second, the total bone shape was obtained by performing a morphological dilation followed by an erosion of the mineralized primary bone segmentation. The perimeter of the total shape is shown as a white line in Video S2. Bone segmentation: while the mineralized embryonic rudiment is not yet truly bone, we use the terminology of bone for brevity. In order to remove noise from the images and improve the precision of the segmentation, a 2D median filter was applied sequentially in all 3 planes using a disk of variable size for each sample as a morphological structuring element adapted to the noise level in each image. Then, two binary volumes of the mineralized parts of the tomographic images were generated (Fig. 1A). First, the mineralized primary bone was selected by contrast threshold (rendered in green in Video S2). Due to higher noise level in few samples (one control TS23, two controls TS24, and three mutants TS24), minor manual adjustments to the segmentation were performed using Dragonfly 3.6, ORS software. Second, the total bone shape was obtained by performing a morphological dilation followed by an erosion of the mineralized primary bone segmentation. The perimeter of the total shape is shown as a white line in Video S2. B) Quantitative analysis: nine parameters were extracted for each sample: Bone volume (BV) and Total volume (TV) values were calculated by summing the pixels in each binary segmented image, from the two segmented volumes (see A) respectively (Fig. 1B1). Maximum bone area (B.Ar max ) and maximum total area (T.Ar max ) were defined as the maximum areas obtained among all the 2D transversal slices, from the binary mineralized primary bone volume and the total bone shape volume respectively ( Fig. 1B2). Among sagittal slices of the total bone shape volume, maximum length (L max ) and minimum length (L min ) were obtained as the longest and shortest distance, respectively, between consecutive white pixels along the z direction ( Fig. 1B3). For each transversal slice, widths along x and y were measured as the longest distances side to side (Fig. 1B2) and eccentricity, calculated as the ratio width x/width y, was used to characterize the cross-sectional bone shape (1 = circular shape, 1 ∕ = elliptical shape). The evolution of T.Ar and B⋅Ar, width x and width y, and eccentricity along the longitudinal length (z axis) were used to visualize and identify the bone tapering. To do so, all images were oriented and sliced along and perpendicularly to the transversal direction. In order to compare among different samples, the maximum length was normalized.

Growth plateB)
C) 3D analysis of chondrocyte properties (Fig. 2D): each chondrocyte was labeled in the volume. The following properties were automatically measured and classified depending on the region to which the corresponding chondrocyte belongs: volume was determined as the sum of the pixels belonging to each chondrocyte times the pixel size. To determine size, the chondrocyte was fitted to an ellipsoid and the length of the 3 axes were measured. The chondrocyte main axes were oriented along to the growth front. Shape: the elongation (ellipsoidity) was calculated as the ratio of the longest axis to the shortest. Orientation was calculated as the elevation angle between the chondrocyte main axis and the horizontal x-y plane. Density was determined as the sum of the chondrocyte pixels divided by the total volume of the zone.
To study the growth plates, a volume of interest was chosen in the center of the growth front of each analyzed rudiment (see Table 1), close to the mineralization front, avoiding the bent regions of the cartilage (see Fig. 2A). The distal growth plates were divided into 4 zones (resting, proliferating, hypertrophic and mineralized zone), based on 2D sagittal slices, as detailed below. Subsequently, a 3D analysis was performed to quantify cellular properties, including cell volume, 3D cell orientation and cell density. Video S3 shows a representative growth plate volume used for the analyses and the image stacks were analyzed using custommade scripts in Matlab.
A) Volume identification and segmentation: the volume of interest (see Fig. 2A) was cropped to select only an 'unbent' region of the growth plate (Figs. 2B1 and S1). The volume was reoriented so that the bone front coincided with the transverse plane (Fig. 2B2). The images were filtered by Contrast Limited Adaptive Histogram Equalization (CLAHE filter) (Fig. 2B3). Contrast thresholding was performed to binarize the images, in which objects were filtered by size to remove noise (Fig. 2B4). As some chondrocytes were touching each other and were jagged in the binary images, one pixel was eroded from the chondrocyte margins in order to separate them (Fig. 2B5). B) 2D analysis of chondrocyte properties to identify the growth plate zones (Fig. 2C) based on orientation, density and shape. Orientation: the chondrocytes were skeletonized to lines and the angle between each line and the x axis was measured. Density: the chondrocyte density was calculated along x (consecutive horizontal lines) as the ratio between white pixels (chondrocytes) and all pixels in line. Shape: the chondrocyte elongation was determined as the ratio between the chondrocytes long and short axes. The results obtained from each z-x slice were averaged over 10 slices and plotted to show the trend along the growth plate length (z axis in Fig. 2) (Fig. S2). The occurrence of each property in the 3D volume was plotted to visualize the distribution of chondrocyte orientation, density and shape along the growth plane length and to identify different zones based on the curve's peaks and valleys (Fig. S2). C) 3D analysis of chondrocyte properties (Fig. 2D): each chondrocyte was labeled in the volume. The following properties were automatically measured and classified depending on the region to which the corresponding chondrocyte belongs: volume was determined as the sum of the pixels belonging to each chondrocyte times the pixel size. To determine size, the chondrocyte was fitted to an ellipsoid and the length of the 3 axes were measured. The chondrocyte main axes were oriented along to the growth front. Shape: the elongation (ellipsoidity) was calculated as the ratio of the longest axis to the shortest. Orientation was calculated as the elevation angle between the chondrocyte main axis and the horizontal x-y plane. Density was determined as the sum of the chondrocyte pixels divided by the total volume of the zone.

Qualitative comparison of phase contrast enhanced tomography and histology
A comparison of chondrocyte and lacunae appearance between two imaging modalities (tomography vs histology) was performed. For histology, forelimbs of two TS24 controls, two TS24 mutants, two TS27 controls, and two TS27 mutants, were embedded in paraffin and longitudinally sectioned at 3 μm thickness using a microtome (Leica Microsystems RM2255, Germany). Humeri were then stained with safranin O. -fast green as described previously [37]. Briefly, tissue sections were deparaffinized using two changes of xylene and rehydrated in ethanol and distilled water. Sections were stained in 0.1% (w/ v) fast green solution and differentiated in 1% (v/v) acetic acid. Sections were then stained in 0.1% (w/v) safranin. Finally, sections were dehydrated in ethanol (95% and 99%) and cleared twice in xylene. Following air-drying, sections were imaged by transmitted illumination using an Olympus BX50 microscope at 4× and 10× magnifications. Chondrocyte sizes were then manually measured using ImageJ. The outcome was compared to the results obtained by the 3D analysis of the chondrocyte properties.

Statistics
A nonparametric Mann-Whitney U test with significance level of 0.05 was performed to compare data from controls and mutants for groups of adequate size. Additionally, a power analysis (G*Power 3.1.9.2 software [38], given α = 0.05, sample size and effect size d) was performed to test the effect of the group sizes.

Bone size and shape are affected by lack of muscular loading during embryonic development
In the case of wild type mice (controls), humerus growth and shaping continued throughout time and by TS24 the bone already displayed an elongated shape similar to the final shape (Fig. 3A). At TS23 and TS24, muscleless limb mice (mutants) had substantially smaller bones with irregular shapes in contrast to the bigger control bones (Fig. 3A, B). In mutants, most of the mineralization occurred after TS24, and by TS27 the size difference between controls and mutants were reduced (Fig. 3A,  B).
Quantitative analysis of the structural changes showed that at TS24, both TV (Fig. 4A) and BV (Fig. 4B) were lower in mutants compared to controls (p < 0.05). Similar tendencies were observed at TS23 and TS27, without being significant. However, the power analysis showed that this could be attributed to the limited sample size in the two groups (Table S1). The same findings were observed for the longitudinal length of the bone (Fig. 4F and G). Similar tendencies were observed for the cross-sectional area without statistical significance ( Fig. 4D and E). Interestingly, the quantitative analysis indicates that at all stages the bone volume fraction (BV/TV) was similar in controls and mutants, this could suggest that lack of muscular contraction primarily affects growth rates but not the relative bone formation (Fig. 4C).
We also investigated how the lack of muscular loading affects bone shape by following the evolution of the total cross-sectional area and the primary bone cross-sectional area along the longitudinal length of the rudiments (Fig. 5). For control mice, the cross-section was wider at the bone's superior and inferior peripheries than in the diaphysis indicating that the bone had a sigmoidal profile. In the case of muscleless-limb mice, bones were more cylindrical, with the typical bone tapering absent at TS23 and TS24, while it was just starting to develop at TS27 ( Fig. 5A and B continuous vs. dashed lines). The difference in width between the epiphysis and diaphysis became more prominent during development ( Fig. 5A and B). The same trend was observed if the moment of inertia was considered, instead of the bone width. The eccentricity values measured for the bone cross-section showed that from TS24, the control bones were elliptical and the ellipticity at the epiphyses increases with further development (Fig. 5C). On the contrary, the cross-section in mutant bones showed eccentricity that remained almost circular along the whole length at all stages.
At TS27 the tapering along the longitudinal axes was visible in control mice (Fig. 6A, yellow line along the muscle profile and Video S4), whereas it was mostly absent in the mutant mice (see Fig. 6A and B, left vs right side, and Video S4). At TS23, no clear differences were apparent in the deltoid tuberosity between controls and mutants since the tuberosity in both controls and mutants was still not distinctly formed.
However, the deltoid tuberosity was largely different between controls and mutants at later stages. In TS24 controls, bone tapering was already visible, the tuberosity had started to form and contained mineralized microstructures (Fig. S3). At TS27, the deltoid tuberosity was well formed in controls whereas in mutants it only started to be visible ( Fig. 6 and Video S4) and its inner volume was largely lacking primary bone (see * in Fig. 6 and Video S4). At TS27 in controls, the deltoid eminence was centered at or above 60% of the normalized longitudinal length while in mutants the eminence was lower, at around or lower than 50% of the normalized longitudinal length (Fig. 6C).

Chondrocyte analysis in controls and mutants at successive embryonic stages
As this represents the first study in which embryonic chondrocytes were imaged by phase contrast enhanced synchrotron tomography, we compared the chondrocytes appearance within their lacunae to reference histology. Our analysis showed that chondrocytes imaged by phase contrast enhanced tomography had similar characteristics (such as size and shape) as when they were stained for histology, whereas the lacunae were less distinguishable in the tomographic volumes (Figs. 7 and S4).
Once that we confirmed that chondrocytes could be clearly detected by phase contrast enhanced tomography, the growth plates of control and mutant embryos were classified into resting zone, proliferative zone, hypertrophic zone and mineralized zone (bone) depending on average chondrocytes 2D size, density and shape. The division into zones was validated on five samples based on visual characterization and comparison by two independent assessors (NN, SA) (Fig. 8A).
The length of the proliferating and pre-mineralized hypertrophic zones was measured (Fig. 8B). In both controls and mutants, the proliferating region appear to increase from TS23 to TS24 and to decrease from TS24 to TS27, while the length of the hypertrophic zone remained quite constant through the embryonic stages. The length of the   6. Deltoid tuberosity and bone tapering are largely absent in muscleless-limb mutants. A) Sagittal slices of phase contrast-enhanced X-ray synchrotron tomograms from a control (left) and a mutant (right). In controls, the deltoid eminence (arrow) is well formed at TS27. Pectoral muscle (P) insertion is evident (arrowhead). The influence of the muscle on tapering along the longitudinal axes is visible (line) along the deltoid muscle profile (D). In mutants the deltoid tuberosity only starts to be visible at TS27 (arrow) and the central volume of the eminence lacks primary bone (*). B) Volume rendering of a control (left) and a mutant (right) rudiment at TS27. C) Segmented total bone area in the normalized longitudinal plane showing the position of the eminence relative to the total bone length (bone lengths were normalized from 0 to 100 so that the relative position of the deltoid tuberosity can be compared between different bones). two zones, and the relative proportion of hypertrophic and proliferating zones, were similar in control and mutants. Importantly, as our aim was to investigate the properties of the pre-mineralized tissue, using phasecontrast-enhanced tomography allowed us to distinguish between premineralized and post-mineralized hypertrophic chondrocytes.
The 3D analysis indicates that in all regions and at all time points, the chondrocyte densities appear higher in controls than mutants, with the greatest difference at TS24 (Fig. 8C). In all regions of the growth plate, the size of chondrocytes decreased with development for both controls and mutants, and mutant chondrocytes seem always smaller than in controls ( Fig. 9A and B). At early stages, the chondrocytes in mutant samples were less elongated (greater short axes/long axes ratio) than control samples in the hypertrophic (at TS24) and proliferating regions (at TS23 and TS24). However, the shape difference disappeared by TS27, when also control chondrocytes displayed lower elongation (Fig. 9C). The qualitative comparison of the outcomes of the semiautomatic analysis (Figs. 8 and 9) to state-of-the art histology showed a good agreement between the techniques (Fig. S4I-K).
Interestingly, unlike the other measured 3D structural properties, the chondrocyte orientation did not significantly change in time or when comparing mutant mice to control mice (Fig. 10). The angles between chondrocytes and the horizontal plane were dispersed in the hypertrophic zone (highlighted in cyan in Fig. 10), reflecting the low orientational organization of the chondrocytes in this area. In the inner proliferating and resting zones, chondrocyte major axis were almost perfectly aligned with the horizontal plane (peak at 0 • , highlighted in yellow in Fig. 10A) while towards the bone edges the chondrocytes tended to follow its curvature and inclination (peak at around ±80 • , highlighted in pink in Fig. 10A). The upper part of the resting zone was partially bent, as reflected by the wider-angle distribution (Fig. 10B). The difference in chondrocyte orientation between the hypertrophic and the proliferating regions, namely from disarranged to similarly-oriented chondrocytes, clearly marked the end of the hypertrophic region.

Discussion
Muscular contractions in utero are essential for normal skeletogenesis [1,2]. Cutting-edge high resolution phase contrast synchrotron X-ray tomography was used in this study to characterize and quantify the effects of lack of embryonic skeletal muscles on developing humeri in 3D, with a micrometer resolution. Differences in bone size and shape, as well as cellular organization in the growth plate were observed between control and muscleless-limb mice.
At early developmental stages, muscleless-limb humeri appeared smaller, with less mineralization, than controls. However, these apparent differences were not statistically significant. Most of the mineralization in muscleless-limb mutants seem to occur after TS24, and by TS27 the size of the mineralized region in mutants was comparable to controls (Fig. 3). This growth trend for controls and mutants is also indicated by the fact that TV, BV, cross-sectional areas and longitudinal lengths appear lower for mutants compared to controls, at earlier stages, with the deficit decreasing with developmental time (Fig. 4). The data also indicate that the main differences between controls and mutants were at TS24. The fact that the size discrepancies between controls and mutants decrease as the fetuses grow, and in particular the recovery of bone size between TS24 and TS27 in the mutant mice, agrees well with hypothesis that mechanical stimuli from the mother and littermates can be transferred passively to the mutant limbs either via flow of the amniotic fluid or through direct pressure [14]. In general, this high resolution 3D evaluation agrees and provides validation for earlier indications based on 2D histology and expand what was previously reported for muscleless-limbs of both Pax3 and for Myf5 − /− :MyoD − /− Fig. 8. Growth plate division in regions and zone lengths. A) The agreement between segmentation based on chondrocyte properties and manual classification was excellent, as exemplified for a control mouse at TS23 (R: resting region, P: proliferating region, H: pre-mineralized hypertrophic region and B: bone region or postmineralized hypertrophic region). B) Lengths of the hypertrophic and proliferating zones for controls and mutants at different development stages. C) Chondrocyte density in the growth plate zones for controls and mutants at each considered embryonic stage (each point represents the mean of all chondrocytes per zone, ranging between 400 and 4000 cells). The lines indicate the mean values and standard deviation. The analysis was conducted for n = 3 controls and n = 3 mutants at TS23, n = 3 controls and n = 2 mutants at TS24, and n = 2 controls and n = 3 mutants at TS27. mice [2,10,11,17,39].
An essential aspect when considering the overall bone shape is the development of the deltoid tuberosity. Bone tuberosity is formed in two steps by a unique pool of progenitor cells that are different from the progenitors that determine the primary cylindrical structure of the bone [40,41]. The first step is regulated by signals from the tendon progenitors and occurs even in absence of muscles. Only the second phase of growth depends on muscular stimuli [40]. Our study expands on what has previously been reported [42,43], by showing that the deltoid eminence is already visible and mineralized at TS24 in controls (Fig. S3), whereas in mutants the tuberosity has not yet formed. It only starts being visible at TS27, but still lacks the internal primary bone microstructure ( Fig. 6 and Video S4). Furthermore, the fact that the tuberosity forms at a lower relative position to the total length in mutants compared to control mice (Fig. 6 bottom) may be linked to the main effect of mechanical loading being in the longitudinal direction of the bone, as indicated by our results (Fig. 4; differences between controls and mutants in F and G vs D and E). These findings are particularly interesting in the light of the work by Stern et al., who showed that in normal conditions long bones grow by isometric scaling [15]. During isometric bone growth superstructures (such as deltoid tuberosity and condyles) maintain their relative positions, which require a constant balance between the proximal and distal growth rates [15]. Even though musclesless-limb mice were not investigated by Eyal et al., their results strongly suggest that muscle and superstructure patterning are regulated in a coordinated manner [44]. Our findings of altered tuberosity position and reduction of longitudinal elongation in mutants could confirm that muscle contractions are playing a fundamental role in keeping a correct growth balance.
As the growth balance seems fundamental for growth and scaling of long bones, and is affected by lack of skeletal muscles, it is highly relevant to establish if muscular loading also affects the growth plate microstructure. Our results show that in control humeri the length of both pre-mineralized hypertrophic and proliferating zones did not differ significantly between mutants and controls (Fig. 8B). This may seem contradictory to the histological results presented by Nowlan et al. showing that the hypertrophic chondrocyte region is reduced in muscleless-limb mice compared to controls [40]. However, in the current study our methodology with phase-contrast-enhanced tomography allowed us to clearly distinguish the microstructures that were already mineralized. We consequently could determine precisely where the mineralization front started and separately consider pre-mineralized and mineralized hypertrophic chondrocytes. This allowed us to determine that even if the bone part is smaller for muscleless-limb mice compared to control mice, the extension of the pre-mineralized hypertrophic zone did not differ. However, even if the lengths of the different regions of the growth plate were similar between controls and mutants, chondrocyte density, size and, in some cases, elongation seem to decrease in response to lack of muscles (Figs. 8C and 9). Interestingly, in the case of controls, the data indicate that chondrocytes reach the highest elongation at TS24, while in mutants the chondrocyte elongation remains more constant at all developmental stages (Fig. 9C). Lower elongation corresponds to flatter chondrocytes which has been suggested to be, at least in part, correlated to the ability of chondrocytes to form columns and consequently to the rate at which the bone is growing [17]. Lower elongation at TS24 could then be a further indication that in Fig. 9. Structural changes at the growth plate during embryonic growth. A) Chondrocyte volumes, B) chondrocyte main axis lengths (parallel to the growth front) and C) chondrocyte elongation (1 = spherical cells) for control mice (blue) and mutants (orange). Each point represents the mean obtained from one sample which was calculated from each region (including 400 to 4000 chondrocytes). The lines indicate the mean values and standard deviation. The analysis was performed for n = 3 controls and n = 3 mutants at TS23, n = 3 controls and n = 2 mutants at TS24, and n = 2 controls and n = 3 mutants at TS27. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) controls, bone growth is fast at this stage. In the mutants, in contrast, most of the growth occurs after TS24, as the biggest difference between control and mutants is at TS24. Size differences are reduced by TS27, possibly due to passive mechanical stimuli compensating for the lack of active movements [14]. Our findings support and expand previous results obtained by 2D manual histological measurements, which showed that lack of muscles leads to a reduction in chondrocyte proliferation [17,45]. However, in our study the properties were measured automatically and in 3D. Thus, our approach has the advantage of reducing bias and errors, and extending the analysis to 3D allows analysis of greater datasets (e.g. for each growth plate we measured the properties of thousands of chondrocytes). For instance, the chondrocyte's 3D orientation was shown to be a good indicator of the passage from the hypertrophic to the proliferating region. It further showed that chondrocytes arrange to follow the bone curvature and inclination both in control and mutant samples (Fig. 10). Moreover, our approach can in the future be extended to analyze chondrocyte column orientations in 3D, to confirm the histological analysis by Shwartz et al. indicating that Splotch delayed muscleless limb mice are characterized by shorter and less aligned chondrocyte columns [17].
Using cutting-edge synchrotron-based techniques restrains the number of samples that can be studied due to limited beamtime, and previous phase contrast synchrotron X-ray tomography studies were typically limited to only two groups with group sizes of 1-4 [22,24,25,27].This study compared 3 time points with group sizes of 4-5, allowing us to even detect some statistically significant differences between controls and mutants in the mineralized rudiment (Fig. 4). The growth plate analysis of cells requires higher image quality (high resolution and complete lack of micro-vibrations during the scan) than the bone analysis, which reduced the number of samples used for the cell analysis to 2 to 3. We conducted both 2D and 3D analyses of the growth plate. The 2D results were used to distinguish different regions of the sample and measure their length, but not to calculate physical chondrocytes properties. For a complete characterization of the chondrocyte properties we performed 3D analysis. The results obtained for the growth plate volumes were compared to manual measurements of histological sections and the outcomes of the two techniques show good correlation (Figs. 7 and S4). Since an extensive visual comparison between the morphological characterization of chondrocytes and lacunae by phase contrast synchrotron tomography and by conventional histology was previously conducted for articular cartilage, for this study we decided to perform the comparison for TS24 and TS27 controls and mutants [46]. However, one interesting point is that lacunae were more distinguishable via histology compared to phase contrast tomography (Fig. 7). In histological images, chondrocytes appear more retracted from the lacunae edge compared to in the phase contrast tomography ( Fig. 7 and D). This could possibly be due to shrinking artefacts during histological processing [47]. This deformation of the chondrocytes has in the past caused many misinterpretations [48]. The fact that the lacunae are less pronounced in the tomographic images could be a further indication that phase contrast enhanced synchrotron tomography has the potential to visualize cartilaginous tissue in a more physiologically representative manner compared to traditional histology. Despite the discussed limitations, one advantage of phase contrast synchrotron X-ray tomography over e.g. using conventionally stained samples and laboratory X-ray tomography is that the staining used to visualize soft tissue by microCT was previously shown to not always efficiently penetrates tissues [20,49]. Moreover, even if the staining penetrates and the absorption of the soft tissue is improved, this often does not lead to improved ability to resolve the internal structures [50]. Consequently, new effective methodologies that allow the study of both mineralized tissue and cartilage are needed. Recently, a step in this direction was presented by Gabner et al., by visualizing both hard and soft embryonic tissues in 3D [51]. The protocol involves imaging Ruthenium Red stained samples by laboratory microCT, which allows a striking visualization of both cartilage and bone at the organ level. However, the bone microstructure and the chondrocytes are not visible. Reaching the resolution required to distinguish microdetails using a laboratory microCT would require extremely long acquisition times (several hours vs. about 10 min when imaging using synchrotron radiation) with Fig. 10. Chondrocyte orientation is not affected by lack of skeletal muscles. A) Angles distribution between the chondrocytes' main axis and the horizontal plane (presented as percentage) B) 3D volume representing schematically the areas of the growth plates that contributed to a certain angle population. consequent possible dehydration of the samples. Furthermore, shrinkage due to dense contrast agents is a commonly reported problem [50].
Our innovative study shows that phase contrast synchrotron X-ray tomography is valuable in the developmental biology field. It gives the possibility of fast 3D visualization of both hard and soft biological tissues at high-resolution (around 1 μm) and prevents misinterpretation of microstructure properties (such as size, shape and orientation) due to the lack of 3D information or due to tissue alteration.
Furthermore, the proposed automated analysis is efficient to characterize large datasets in shorter time and to reduce bias. The workflow could be adapted to study other complex biological tissues e.g. bonetendon interfaces; periodontal ligament, dentin and enamel in teeth; mineralization in shells and planktonic animals; silica, calcium and lignin deposits in plant tissues.

Conclusions
High-resolution phase contrast X-ray tomography was used to increase our understanding of the mechanobiology of bone formation in embryonic long bones and its dependence on muscular loading. The proposed protocol can capture the 3D morphology of both mineralized and soft tissues. We demonstrated that, for the mineralized region and the growth plate, the absence of muscular contractions has particularly pronounced effects at TS24. The impact of lack of muscular loading on bone shape development at later stages (TS27) is characterized by more cylindrical bones, impaired elongation and altered relative position of the tuberosity. In the primary growth plate, if the lengths of the different regions did not seem impacted, the absence of muscle contractions seem to result in lower chondrocyte densities, with cells of reduced sizes and less elongated.