Multiscale analysis of 3D nuclear morphology reveals new insights into growth plate organization in mice

The shape of the nucleus is tightly associated with cell morphology, the mechanical environment, and differentiation and transcriptional states. Yet, imaging of nuclei in three dimensions while preserving the spatial context of the tissue has been highly challenging. Here, using the embryonic tibial growth plate as a model for cell differentiation, we study nuclear morphology by imaging cleared samples by light-sheet fluorescence microscopy. Next, we quickly segmented tens of thousands of nuclei using several open-source tools including machine learning. Finally, segmented nuclei underwent morphometric analysis and 3D spatial reconstruction using newly designed algorithms. Our method revealed differences in nuclear morphology between chondrocytes at different differentiation stages. Additionally, we identified different morphological patterns in opposing growth plates, such as gradients of volume and surface area, as well as features characteristic of specific growth plate zones, such as sphericity and orientation. Altogether, this work supports a link between nuclear morphology and cell differentiation. Moreover, it demonstrates the suitability of our approach for studying the relationships between nuclear morphology and organ development. Author summary There has been a growing interest in the relationship between nuclear morphology and its regulation of gene expression. However, to study global patterns of nuclear morphology within a tissue we must address the problem of acquiring and analyzing multiscale data, ranging from the tissue level through to subcellular resolution. We have established a new pipeline that enables acquisition and segmentation of hundreds of thousands of nuclei at a resolution that allows quantitative analysis. Moreover we have developed new algorithms that allow superimposing morphological aspects of hundreds of thousands of nuclei onto a visual representation of the entire tissue, allowing us to study nuclear morphology at an organ level. Using mouse growth plates as a model for the relationship between nuclear morphology and tissue differentiation, we show that nuclei change different aspects of their morphology during chondrocyte differentiation. Growth plates are usually described generically in the literature, suggesting they lack unique characteristics. We challenge this dogma by showing that morphological features such as volume distribute differently in opposing growth plates. Altogether, this work highlights the possible role of nuclear shape in the regulation of cell differentiation and demonstrates that our approach enables the study of nuclear morphology patterns within a tissue.


Introduction
The ability to analyze multiscale data is essential for investigating the contribution of cellular components to tissue morphogenesis. Analyzing three-dimensional (3D) data becomes complicated when they span multiple scales, as is the case of the spatial relationships between cells and tissue. Focusing on the tissue level results in loss of subcellular resolution, while focusing on smaller regions at high resolution results in loss of their relationship to 3D tissue morphology. To cover all these scales, the entire dataset must be at the highest desired resolution, which produces even more challenges regarding data size. To study the cellular contribution to tissue morphogenesis, we need a way to easily explore large datasets describing the tissue level, while preserving individual cell features.
To image whole tissues and visualize in 3D subcellular tissue components, optical sectioning is often applied. Using fluorescence microscopes such as confocal or two-photon, this method produces high spatial-resolution datasets, but at the cost of low imaging speed. Moreover, these techniques are limited in imaging depth and are most suitable for tissues up to 300 µm in thickness. Light-sheet fluorescence microscopy (LSFM), which quickly acquires large volumes of tissue at submicron scale, is ideal for imaging fluorescent tissues or even whole organisms such as zebrafish (1,2), C.elegans (3,4), or Drosophila (5). Whole organism imaging with LSFM is possible because these model organisms are transparent and small, a critical factor for this modality. Thus, applying LSFM to opaque mammalian tissues requires tissue clearing to enable in-depth imaging.
Biological samples absorb and scatter light. Tissue clearing is the process of reducing light scattering and absorption by making tissues optically transparent. This process, which typically involves removing macromolecules and matching refractive indices within samples, results in better penetration of light through the sample, thereby increasing imaging depths (6). Recently, several organs from adult mice were successfully cleared using aqueous-based buffers, allowing for imaging up to 1.5 mm in depth while maintaining endogenous fluorescence (7)(8)(9). Combining tissue clearing with LSFM will allow for high-resolution mapping of fine tissue architectures in an intact 3D sample (10).
For many years, researchers have used information about nuclei to extrapolate about the state of a tissue. Changes in nuclear morphology have been associated with malignant transformation of cells and tissue (11,12), aging (13,14), cell differentiation (15)(16)(17)(18), and response to mechanical forces (19)(20)(21). In the last 15 years, accumulating data have connected nuclear morphology with molecular mechanisms that involve gene expression or transcription (19,(22)(23)(24). Chromatin is organized into open and closed domains, which allow for the regulation of gene expression (25).
The finding that chromatin associates with nuclear lamins and with the inner nuclear membrane provides a mechanism whereby alterations in nuclear morphology can regulate gene expression (26). These findings highlight the importance of studying global patterns of nuclear morphology within a tissue.
The growth plate is an excellent model system for studying the multiscale relationships among nuclear morphology, cell differentiation and tissue geometry. Growth plates are cartilaginous tissues located at either end of developing bones and responsible for their elongation. Starting in embryonic development, this tissue comprises chondrocytes that display a conserved spatial organization into zones, namely the resting (RZ), proliferative (PZ), prehypertrophic (PHZ), and hypertrophic zones (HZ). These zones, reflecting a series of differentiation states, are marked by unique cell morphologies, extracellular matrix (ECM) properties and gene expression profiles (27)(28)(29)(30)(31).
Despite these advantages, composing a multiscale dataset from growth plates poses several challenges. Growth plates are extremely dense, as they are packed with collagen-rich extracellular matrix (31), which hinders 3D imaging. Additionally, the tissue is relatively large, reaching >1 mm in thickness, which implies longer durations of imaging and data analysis.
Finally, current methodologies to segment nuclei are not suitable for large datasets, within which image quality may vary, resulting in slow and unreliable segmentations. In this work, we address these challenges using an especially designed pipeline for 3D imaging and analysis of growth plate cell nuclei as well as their spatial relationships within the tissue. We imaged nuclei from cleared embryonic growth plates using LSFM, applied several open-source tools to quickly segment nuclei, and developed a visualization platform to explore our obtained data. Our method revealed that during chondrocyte differentiation, nuclei underwent several morphological changes that were common to all growth plates, such as in nuclear orientation. By contrast, we identified different patterns in the two opposing growth plates in features such as nuclear volume and surface area.

Results
A pipeline for imaging the growth plate and quantifying large multiscale data To study nuclear morphology in an intact mouse growth plate, we developed a high-throughput method that allowed us to characterize nuclear morphology quantitatively while maintaining the spatial relation of each nucleus within the tissue (Fig 1). To overcome the challenge of imaging a large and dense tissue like the growth plate, we cleared the tissue using the simple and inexpensive PACT-deCAL technique (7,8), which was found to allow nuclear labeling while preserving tissue fluorescence. Next, embryonic day (E) 16.5 DAPI-labeled tibias were imaged by light-sheet microscopy (Fig 2A). Due to its high speed, light-sheet microscopy proved to be the superior choice over other laser-based techniques, such as confocal and two-photon microscopy. To ensure that light would pass through the sample, we used dual side illumination when necessary and rotated the sample to find the optimal imaging orientation. To reduce the file size and expedite image processing, acquired z stacks included only regions with growth plate nuclei.
The next challenge was to conduct accurate and time-efficient 3D segmentation of nuclei from the light-sheet images. For that, we combined three open-source platforms to perform semiautomatic segmentations (SAS). Nuclei from surrounding tissues were manually excluded by masking using Microview 2.1.2 (GE Healthcare). Next, the cleaned raw images were automatically segmented using XPIWIT software (32,33). Then, an object classifier to identify and remove incorrect segmentations was trained using Ilastik (34). 30 segmented images from various regions within the growth plate were used in the training dataset. We were then able to classify hundreds of remaining images blindly by batch processing, thereby greatly speeding up the cleaning stage ( Fig 2B).
One of the consequences of imaging a whole tissue was that the border of nuclei residing in the center of the tissue could appear fuzzy due to lower signal-to-noise ratios (Fig S1.), which could obviously influence the accuracy of subsequent analysis. To overcome this potential problem, we performed a benchmark analysis comparing segmentations from two different imaging modalities, namely confocal and light-sheet microscopy. To establish the "ground truth" (GT), nuclei from the four zones of the growth plate, namely the resting (RZ), proliferating (PZ), prehypertrophic (PHZ), and hypertrophic zone (HZ), were manually segmented from confocal images in Microview. Volume and surface area were extracted and corresponding histograms were compared between manual and semi-automatic segmentations using chi-square distance.
Results showed that nuclear segmentations by SAS were indistinguishable from manual GT segmentations (p > 0.05), demonstrating the validity of our measurements (Fig 3A). Since the segmentation was performed to the same quality on four different growth plates from two separate litters (Fig 3B), we could use the segmentation error as a tool to estimate the total number of growth plate nuclei at this developmental stage (see Methods). We estimated the number of nuclei to be 178,825±35,582 in the proximal growth plate and 90,460±16,448 in the distal growth plate (Fig 3C).

Morphometric and spatial analysis of growth plate nuclei
Having segmented successfully nuclei from each growth plate, we proceeded to conduct morphometric analysis on these nuclei, which included calculating the average volume (Fig 4A), surface area (Fig 4B), sphericity (Fig 4C), and spatial orientation (Fig 4F). To understand the spatial distribution between different nuclei, we calculated nuclear density ( Fig 4D) and occupation ( Fig 4E). Finally, to explore the relationship between nuclear morphology and its location in the growth plate, we stitched the data together and generated 3D morphology maps.
For that, nuclei were mapped back to their original locations in the growth plate by using the centroid of each segmented nucleus and the stage coordinates of each image. To analyze data pertaining to tens of thousands of nuclei, we reduced the data dimensionality by splitting the 3D images into a regular grid composed of 75-µm 3 cubes.
In the proximal growth plate, nuclei gradually increased their volume and surface area up to 9fold as they differentiated. Interestingly, however, nuclei in the distal growth plate did not form a clear gradient. Instead, there were two regions, in between the presumptive RZ and PZ and in the HZ, where nuclei displayed increased volume and surface area compared to the rest of the tissue.
Moreover, maximum volume and surface area were on average 17% and 10% less than in the proximal side, respectively. Relying on previously published data on chondrocyte volume during differentiation (35,36), we found that nuclear to cell (N/C) volume ratios increased from 0.18 in the RZ to 0.39 in the HZ, suggesting that these ratios are not constant during chondrocyte differentiation in the growth plate. Nuclei located on the edges of the growth plate along the medial-lateral axis consistently had lower volumes and surface areas. Previous studies showed that cells in the RZ and HZ exhibit equally high sphericities (35,36). By contrast, the highest nuclear sphericities were found in the PHZ and HZ ( Fig 4C).
Next, we evaluated nuclear density in the growth plate ( Fig 4D). In the proximal growth plate there was a constant gradient along the proximodistal (P-D) axis with highest densities in the RZ and lowest densities in the HZ and along the edges of the growth plate. In the distal growth plate the gradient had two separate phases, a constant increase and a constant decrease along the P-D axis. To study the ratio between nuclear volume and total tissue volume, we calculated the nuclear occupation percentage (Fig 4E). In correlation with our density findings, areas with high nuclear density such as in the RZ had high nuclear occupation and vice versa. To investigate whether different growth plate zones are characterized by unique nuclear morphology features, we performed automated classification of individual nuclei (see Methods) ( Fig 5A).
Interestingly, the spatial distribution of the classified zones recapitulated the normal organization in the growth plate with 74% accuracy, suggesting that the information encoded in nuclear morphology is sufficient to define a growth plate zone ( Fig 5B).
Previous studies showed that cells in the proliferative zone orient themselves into columns that are parallel to the long axis of the bone, facilitating bone elongation (37)(38)(39)(40). Using our segmented data, we directly quantified the orientation of the three major axes of each nucleus using principle component analysis (PCA), allowing us to visualize the spatial distribution of nuclear orientations in the growth plate ( Fig 4F). In accordance with the current knowledge of cell orientation, we found that nuclei in the PZ also oriented their shortest axis (PC3) along the long axis of the bone with little variability. This trait was conserved in other regions of the growth plate, yet with higher variability (Fig S2). By contrast, a high degree of variability in nuclear orientation direction was observed in the HZ and some areas in the RZ. Finally, nuclei located on the edge of the growth plate oriented orthogonally to the surface of the growth plate.

Discussion
Data that describes tissue architecture using only one scale, namely tissue or cells, provides limited information. Without the ability to visualize subcellular morphology and location, it is impossible to identify resulting pattern formation within a tissue. This multiscale approach has recently been utilized in studying epithelial tissue morphogenesis in Drosophila wings (41) and spheroid spatial heterogeneity (42). In these works, multiscale analysis revealed a role for convergent extension in shaping wing veins, and showed that in breast carcinoma spheroids, there are differences in cell densities between different layers, which should be considered and further explored when modeling this type of cancer. Similar to these recent studies, the method we describe here was developed to gain a system-level understanding of tissue architecture through multiscale analyses of nuclear morphology in the growth plate. By combining tissue clearing and LSFM with high-throughput nuclear segmentation, we were able to implement our newly developed algorithms to create high-resolution 3D tissue maps and extract patterns of nuclear morphology and spatial orientation within the large and complex embryonic growth plate.
For high-throughput extraction of nuclear morphology, our data needed to be acquired at high resolution. Imaging at high resolution creates technical challenges, such as longer imaging time and larger resulting datasets. By acquiring z-stacks consisting only of relevant ROIs, imaging speed was increased and data size was reduced. Even so, our datasets were 500-800 GB in size per growth plate. To overcome this problem, we performed down-sampling of each image. As the voxel dimensions in the raw image (0.1 µm) are ~60 times smaller than the average nucleus diameter (6 µm), isotropic down-sampling by a factor of 2 can be done with negligible loss of morphological accuracy. Several rounds of down-sampling performed before and after segmentation resulted in 80% reduction in data size, greatly speeding up the subsequent preprocessing and segmentation steps.
In recent years, many segmentation tools have become openly available. To overcome the problem of speed and accuracy, we performed automatic segmentation in XPIWIT followed by automatic classification of correct and incorrect segmentations using a trained classifier in Ilastik to quickly extract the 3D morphology of tens of thousands of nuclei per dataset. Although manually training the classifier initially took between one and two weeks for one growth plate, this same training was applied on other untrained growth plates, resulting in classification of the remaining growth plates in one day. As the entire segmentation process is semi-automatic, the user still needs to do some work, such as removing regions in the image that do not belong to the growth plate and tuning segmentation parameters if necessary after evaluating the segmentation output. Additional speed and automation could be gained by initiating these processes directly through computer commands.
Accurate segmentation of thousands of growth plate nuclei allowed us to investigate the relationship between patterns of nuclear morphology and the growth plate. It has been shown in previous studies that there is a correlation between cell and nuclear morphology (43)(44)(45). In Drosophila egg chambers, nuclear volume strongly correlates with cell volume (46).
Additionally, in growth plate explants, shape factors such as width/height ratio as well as volume are correlated between cells and nuclei (47). Previous studies showed that chondrocyte volume increases by up to 6-fold while they move from the resting zone to the hypertrophic zone (35,36). In agreement with these studies, nuclear volume increased as the cells moved along the growth plate; yet, surprisingly, nuclear volume increased up to 9-fold. Numerous works have shown that nucleus size is tightly linked to cell size. In both growing yeast cells and Xenopus laevis, the nucleus to cell (N/C) ratio is constant across cells differing in size (48)(49)(50)(51).
Interestingly, our measurements for 50% of all growth plate nuclei revealed distinct N/C ratios in different zones. This finding may suggest that in the growth plate the mechanism that preserves the N/C ratio is variable between the different zones, or that it is not functional during chondrocyte differentiation.
Looking at the orientation maps, we were able to identify the presumptive proliferative zone in both growth plates, where cells stereotypically form columns parallel to the long axis of the growth plate (38). In fact, we found that most nuclei in all zones aligned their shortest axis parallel to the main axis of the growth plate, although the degree of variability was noticeably higher in regions other than the proliferative zone. These findings are in line with previous studies of chondrocyte orientation where it was shown that cells in the proliferative zone orient their longest axis perpendicular to the long axis of the growth plate (52)(53)(54).
Although the two growth plates of the same bone appear to have similar structures, we could identify some morphological variation between them. For instance, while similar patterns of nuclear density, occupation, and orientation were observed, there were differences in volume, surface area, and sphericity. While there was a constant gradient of volume and surface area through the center of the proximal growth plate, the distal end displayed a gradient with a constant increase and constant decrease along with reduced volume and surface area relative to the proximal side. Past studies have outlined morphological and functional differences between the proximal and distal growth plates of different long bones, particularly in overall size and shape and in growth rate (55)(56)(57)(58)(59). The different growth rates between proximal and distal growth plates were recently shown to serve a functional purpose, allowing isometric scaling of symmetry-breaking elements during long bone development (60). Additionally, ATAC-seq on proximal and distal femoral growth plates revealed differences in chromatin accessibility (61).
These findings support the notion that in addition to the generic mechanisms that regulate all growth plates, there are also unique features that are regulated separately in each growth plate.
As nuclear morphology is tightly linked to chromatin accessibility (62) and gene regulation (63), it is plausible that the changes in nuclear morphology observed in different areas of the growth plate reflect changes in gene expression. Thus, our method may provide a new avenue to study the relationship between tissue structure and gene expression in any tissue through the use of 3D mapping of nuclear morphology.

Light-sheet microscopy
The cleared samples were imaged using a light sheet Z1 microscope (Zeiss Ltd.

Image segmentation
Prior to segmentation, images were down-sampled in the X,Y direction using the Downsample plugin in Fiji (64) resulting in x,y,z voxel size of 0.137 µm, 0.137 µm, 0.387 µm. Noncartilaginous nuclei were manually removed using Microview 2.1.2 (GE Healthcare). Shortly, contours were drawn around non-cartilaginous nuclei, and the resulting ROI was reassigned a pixel intensity value of 0. To automatically segment fluorescently-labeled nuclei in the 3D images, we performed a two-step procedure. In the first step, seed points that were roughly located in the center of the nuclei were detected using a Laplacian-of-Gaussian-based (LoG) approach as described in (33). In brief, the 3D input images were filtered with a LoG filter using a standard deviation that was empirically tuned to the radius of the objects of interest. We used Fiji and images were saved as a 3D tiff stack.

Nucleus feature extraction
Following image segmentation, we performed an initial data cleaning step that included removal of nuclei that overlapped with the border of the imaged section, as well as nuclei with irregularly small (< 100 µm 3 ) or large (> 4000 µm 3 ) volumes. To extract morphological features, each nucleus was first converted from a binary volume into a triangle mesh by applying Gaussian smoothing filter (sd = 0.5 pixels) and extracting the iso-surface at the iso-value 0.5. From the triangle mesh of each nucleus, the following features were extracted: surface area, as the sum of the areas of all mesh faces; volume, using the convergence theorem (a.k.a. Gauss's theorem or Ostrogradsky's theorem (65)); sphericity, calculated as: where is the volume of the nucleus and is the surface area (66)

Machine learning classification
To evaluate how well nuclear morphology could predict to which zone a cell belongs, we developed an automated classification method using MATLAB R2018a's "Classification Learner toolbox". We used classification models on manually labeled nuclei based on their volume, surface area, sphericity and density to train a series of classifiers (Table S1). The feature values were first normalized as a score function (i.e., subtraction of the mean value and division by the standard deviation obtained from the entire set of values). The accuracy of each of the classification models was assessed using a 5-fold cross-validation method. The best model for the classification of the nuclei in the proximal growth plate was a Fine Gaussian SVM, with an accuracy of 74%, whereas the best model for the classification of the nuclei in the distal growth plate was Boosted Tree with an accuracy of 69.5%. Confusion matrices describe the details of the prediction. Using the learned classifiers, we made predictions on the unlabeled nuclei.

Segmentation Validation
To validate segmentation, an E16.5 tibia from a sample littermate was dissected and cleared using PACT-deCAL, stained with DAPI, mounted in RIMS on a glass slide with a shallow well, and covered with a glass coverslip. 75-µm z-stacks from areas in the RZ, PZ, PHZ, and HZ were acquired with x,y,z voxel sizes of 0.1 µm, 0.1 µm, 0.5 µm. Confocal images were manually segmented (n = 80 nuclei) in Microview software and volume and surface area were calculated for each nucleus. Volume and surface area of automatically segmented nuclei (n = 80) from comparable regions in the growth plate were compared to the manual segmentations by a chisquare distance between histograms test (67). We could not reject the hypothesis that the two distributions are drawn from the same underlying distribution (p > 0.05), demonstrating that our semi-automatic segmentations were valid.

Segmentation error calculation
To calculate the segmentation error from six representative regions in the proximal and five representative regions in the distal growth plate, we calculated the average ratio of analyzed nuclei out of the total nuclei number in each region of 2D optical sections every 23 µm along the Z axis, to avoid counting the same nucleus twice. To count the total nuclei number in each slice, we counted nuclei in the DAPI channel, using a macro written in FIJI. To count the number of analyzed nuclei out of each slice, we used the MATLAB scripts. The resulting calculations were done using the following equations: Average error = total nuclei / analyzed nuclei Estimated total nuclei = # analyzed nuclei * (Average error)

STD = std (Average error)
Average segmentation error = 1 / average error per growth plate

3D nuclear morphology maps of growth plates
To visualize the spatial distribution of nuclear morphology within each growth plate, data were first stitched back together in MATLAB using the stage coordinates and scaling of each image.
To place each voxel in a global coordinate system, we added the coordinates of each stack to every voxel in the stack, thereby reconstructing the entire bone. To highlight large patterns while averaging out small differences between individual nuclei, we chose to represent nuclear features at a coarse-grained level. This representation was computed by first defining a regular spatial grid over the data. Each element of the grid was defined as a 75 x 75 x 75 µm 3 cube, containing up to 275 nuclei. Within each cube, we averaged the nucleus volume, surface area and sphericity and used the Delaunay tessellation field estimator at a resolution of 2 µm 3 to compute the nuclear density for each cube (68). In addition, for each cube, we computed the occupation, defined as the sum of nuclear volumes divided by the volume of the cube. Then, using MATLAB's jet color map, we represented the characteristics of each cube on the grid by drawing spheres whose radii and colors are proportional to the computed values and whose centers correspond to the average position of the nuclei within the bin.

Nucleus orientation maps and scatter plots
The spatial distribution of nuclear orientations was visualized using a regular spatial grid of the To compare among growth plates, we registered the stitched data beforehand. The registration was performed in two steps using MATLAB's rigid registration. First, each growth plate was registered to an artificially generated set of points mimicking the global shape of a growth plate, with the P-D axis aligned with the X axis. Then, each growth plate was registered to one of the samples. For each cube a line was drawn oriented and color-coded according to the average orientation, and whose thickness was proportional to the spherical variance. Finally, to summarize the distribution of orientations in a 2D graph, we used the fact that each of the mean orientation vectors obtained for each of the bin of the regular grid were encoded using spherical coordinates ( and , for polar and azimuthal angles respectively). varies from 0 to , while varies from − to , enabling the encoding of any orientation in 3D.

Spatial profiles of morphology features
To visualize quantitatively the spatial profile of each of the computed features and identify gradient along the growth plate, we defined an irregular grid as following: The entire growth plate was divided into 75 µm-thick slices along the P-D axis, and into three equally spaced slices along each of the two remaining axis. This division resulted in 9 rods organized in a 3      Figure S1. Signal-to-noise ratio is reduced away from the growth plate surface (related to Figure 3). Optical sections from LSFM imaging of cleared growth plates with their accompanying intensity histograms generated in Fiji (below) demonstrate that the signal-to-noise ratio in nuclear images is lower in the center of the tissue than on the edge. Figure S2. Spatial distribution of nuclear orientation (related to Figure 4). Scatter plots of nuclear orientation defined by polar angle and azimuthal angle for the proximal and distal growth plates. Colored dots represent distribution of the colored lines in Figure 4F. Large concentration of green dots highlights the existence of a preferential nuclear orientation in both growth plates. Table S1. Proportions of manually labeled nuclei in each growth plate zone (related to Figure 5). Automatically segmented nuclei from the proximal and distal growth plates were labeled as belonging to the HZ, PHZ, PZ, or RZ. The volume, surface area, sphericity, and density of these labeled nuclei were used to train a series of classifiers to predict which zone they belong to.