Automated segmentation of thick confocal microscopy 3D images for the measurement of white matter volumes in zebrafish brains


 Tissue clearing methods have boosted the microscopic observations of thick samples such as whole-mount mouse or zebrafish. Even with the best tissue clearing methods, specimens are not completely transparent and light attenuation increases with depth, reducing signal output and signal-to-noise ratio. In addition, since tissue clearing and microscopic acquisition techniques have become faster, automated image analysis is now an issue. In this context, mounting specimens at large scale often leads to imperfectly aligned or oriented samples, which makes relying on predefined, sample-independent parameters to correct signal attenuation impossible.
 Here, we propose a sample-dependent method for contrast correction. It relies on segmenting the sample, and estimating sample depth isosurfaces that serve as reference for the correction. We segment the brain white matter of zebrafish larvae. We show that this correction allows a better stitching of opposite sides of each larva, in order to image the entire larva with a high signal-to-noise ratio throughout. We also show that our proposed contrast correction method makes it possible to better recognize the deep structures of the brain by comparing manual vs. automated segmentations. This is expected to improve image observations and analyses in high-content methods where signal loss in the samples is significant.


Introduction
Due to its small size, its robustness, affordability and its high reproduction rate, zebrafish have become a popular vertebrate model over the last decades, which is used in a variety of screening experiments [7]. These screens are performed for the study of behaviour [4], gene expressions [10], drug-toxicity [16] or exposure to toxicological compounds [3]. However, image-based screening approaches, which are using whole zebrafish embryos, are mainly performed on the basis of 2D images [17,19]. While 3D-based screening methods do exist, they are either using unconventional microscopy setups, or need multiple acquisitions of the same sample [2,8,14,15,21], which renders them unsuitable for large scale analysis.
Recently developed tissue clearing methods allow biologists to perform anatomic studies on the basis of in-toto confocal image acquisitions of thick samples, such as mouse brain or entire zebrafish [6]. However, tissue clearing of large samples is time consuming. In addition, the most effective techniques involve hazardous compounds like organic solvents [20] which complicates their application at large scale. Furthermore, some tissue clearing protocols have the potential of changing the size or -even worse -the shape of the processed tissue [23]. Simpler protocols are now becoming available; they are reducing size alteration, and are both less toxic and faster [1].
These fast tissue clearing protocols open the route to new high content screenings (HCS) using deep confocal imaging on tissue cleared samples, like whole-mount zebrafish larvae. Such protocols allow the detection of subtle abnormalities of larvae due to drug exposure or mutations. This approach is based on two technological developments: first, fast mounting procedures, which are an improvement of earlier work [22] and allow the reliable mounting and imaging of numerous samples. Second, the automatic image analysis, combining signal improvement, feature detection and measurement of regions of interest.
In this article, we describe an automatic and accurate segmentation procedure for measuring the volumes of whole larvae and the volumes of brain white matter of in-toto 5 days post-fertilization (dpf) zebrafish imaged by confocal microscopy. Our protocol involves a new and fast contrast correction step, which reduces the impact of light attenuation and dispersion in the samples.
Indeed, current contrast corrections are based on the Beer-Lambert attenuation model [13] which demands time-consuming computations and often require multiple acquisitions in different orientations to achieve robust results [12,18]. Using the median value at each depth of the acquired sample, we show that our method reproduces a Beer-Lambert attenuation over the sample while accelerating computation by reducing the algorithmic complexity of the contrast correction.

Tissue clearing
Zebrafish larvae are fixed in formaldehyde (4% formaldehyde, 0.1% Tween20 in 1X phosphate buffered saline (PBS)) and washed in 0.1% Tween20/PBS (PBStw). Fixed larvae are then bleached by incubation in H2O2 -based depigmentation solution for 2.5 hours (1.5 hours in 0.5X saline sodium citrate buffer (SSC), 0.1% Tween20 in dH2O followed by 1 hour in 0.5X SSC, 5% formamide and 3%H2O2 in dH2O) and then thoroughly washed in PBStw. Larvae are dehydrated in a series of ethanol concentrations (25%, 50%, 75% and 100% in PBStw) and stored at -20°C in 100% methanol. After rehydration in a series of ethanol concentrations (100%, 75%, 50% and 25% in 0.5% Triton X-100 in 1XPBS (PBStr)), larvae are washed in PBStr prior to proteinase K treatment (at 20 µg/ml in 0.5% PBStr). They are post-fixed (2% formaldehyde, 0.1% Tween20, DMSO 2% in 1X PBS) for 1 hour and washed in PBStr. Larvae are blocked 4 hours in 10% NGS, 0.05% Azide in 1XPBStr (IB) and finally labeled with 1 mM Dil (DiIC18(3) Stain, Molecular Probes) in IB for 2 days and washed in PBStw. We used the fluorescent lipophilic DiI stain, which causes an intense and consistent staining predominantly in the white matter due to the high density of lipid-rich myelinated fibers. Due to the comparatively small size of the specimens, the applied tissue clearing procedure only involves the matching of the refractive index (RI) between the specimen and the surrounding medium. In the context of HCS on 5 dpf zebrafish larvae, this protocol is currently the fastest, safest and most economic way for generating a large number of transparent specimens. The tissue clearing medium is a fructose-based high-refractive-index solution (fHRI) prepared as follows: 70% fructose (wt/vol), 20% DMSO (wt/vol) in 0.002 M PBS, 0.005% sodium azide (wt/vol). The RI of the solution was adjusted to 1.457 using a refractometer (Kruss, Germany). Two kinds of mountings are used. To evaluate accuracy of white matter segmentation, specimens were mounted with dorsal side up in a Petri dish by embedding them in a drop of 0.8% Phytagel (Sigma-Aldrich, #P8169) in embryo medium. To evaluate signal decrease over depth of the sample, in particular below the eyes, specimens were mounted between two coverslips in lateral position to allow imaging from both the left and right sides. To prevent motion between both side acquisitions, specimen were embedded in 0.8% low-melting agarose within 700 µm thick silicone spacers. The RI is equilibrated by incubating the specimen in fHRI for at least 12 h and another coverslip is then added and sealed with silicone.

Confocal microscopy
Confocal microscopy image stacks (3D images) of whole-mount zebrafish larvae (5 dpf) were acquired with a 12-bit encoded dynamic range. DiI signals were recorded using a Leica TCS SP8 laser scanning confocal microscope equipped with a Leica HC FLUOTAR L 25x/1.00 IMM motCorr objective. For specimens embedded between two coverslips, acquisitions of opposite orientations of the specimens are registered by a two-step process using the Elastix software package [11]. First, the so-called moving image is flipped horizontally and the order of its slices is reversed to compensate for the physical flipping of the specimen. Second, the moving image is aligned to the fixed image in a multi-resolution fashion (using 4 resolution levels of a Gaussian pyramid in our case), using an affine transformation to model the deformation between both images, and either the Mean Square Difference (MSD) or Mutual Information (MI) of the voxel intensities is used as similarity metric. Finally, the transformed moving image and the fixed image are fused using their maximum voxel intensities.

Segmentation accuracy
We manually performed 3D segmentations of whole larvae and their white matter using using 3D Slicer [9] and Amira for Life & Biomedical Sciences (Thermo Fisher Scientific) software packages. These manual segmentation were then compared to corresponding automatic ones.
We measure the segmentation accuracy with a Matthews correlation coefficient (MCC), which considers both true positives and true negatives. We also use Sørensen-Dice coefficient (DSC), which does not consider true negatives. To balance the importance of false positive over false negative voxels, we also used the general balanced metric [5]. Setting δ = 2, this metric is called GB2 in the sequel. Coupled with DSC, this metric informs on the ratio of false positive voxels over false negatives ones.
The imaging procedure leads to numerous black pixels, which could artificially impact positively the MCC. To avoid this artefact, we chose to reduce each acquisition to the bounding box of the union between the ground truth and the result of the segmentation, with a tolerance of 5 voxels in each direction. This bounding box is computed to reduce the amount of true negative voxels, thus reducing their impact on the MCC computation.

Algorithms
The following sections use notations given in Table 1.

Whole larva segmentation
Let I i be the initial 3D acquired volume.

Notations
Name Definitions To initiate whole larva segmentation, a morphological opening using a cubic structuring element of size 15 × 15 × 15 voxels ¹ is performed to smooth the signal (Eq. (1)). Then, a first mask of the segmentation is computed. To achieve this, the opened image is thresholded (Eq. (2)) using the first percentile of grey values of non-black voxels of I i (Eq. (3)) as thresholding value. Then, a marker image is created by setting local maxima of I i inside A 2 (Eq. (4)) to 2, the frame of the image to 1 and every other voxels to 0 (Eq. (5)). A morphological gradient is then performed on the opened image (Eq. (6)). Then, a classical watershed is computed using the previously explained marker image and the morphological gradient (Eq. (7)). This resulting image is closed using a cubic structuring element of size 11 (Eq. (8)) to smooth the resulting volume. The last step of the segmentation of the whole larva is an ultimate volume opening (Eq. (9)). This step consists of performing a succession of volume openings until a single object remains on the array. This last object is kept as the final larva segmentation result.

Segmentation-based depth-dependent contrast correction
Using the previously obtained segmentation of the whole larva S larva , we now compute the depth map ( Fig. 1B) of the specimen along the Z axis. The first step is a cumulative sum of S larva along the Z axis (Eq. (10)). Then, we compute the multiplication of the original segmentation and the cumulative sum, to reduce values to the previously segmented shape, as in Eq. (11). The results of this multiplication is the depth of each pixel into the specimen starting from the acquisition side.
Let d be the value of the depth in the sample, d b the last value of d, µ(d) the median grey value of the original image at level d, µmax the maximum value of µ(d) and dmax the index of µmax in µ(d). For each d from dmax to d b we compute the new voxel value by multiplying each grey value of d by µmax divided by µ(d) (Eq. (12)).
A result of this procedure is displayed on Fig.1

White matter segmentation
The first step of the segmentation of the grey matter in our sample is the segmentation-based depthdependent contrast correction computation of I DiI . Let C DiI be the contrast corrected image. We next compute a median filter with a cubic structuring element of size 3 on C DiI (Eq. (13)). After this median filter, we compute a four-class pixel binning based on Otsu thresholding, from which we retain only the two highest classes (Eq. (14)). This step provided a rough segmentation of our region of interest. However, this segmentation also contained eyes and yolk. We now propose additional steps to remove these tissues.
Due to the DiI labeling, the retina appear over-saturated on the image. Therefore, we use this information to detect retinas with a threshold 99% of 2 12 − 1(Eq. (15)). This value was chosen because our datasets were  16)). Taking every pixel brighter than 99% of 2 12 flags saturated pixels and their neighborhoods. A volume opening of size 200 is then computed (Eq. (17)), and the result of this opening is closed using a cubic structuring element of size 7 (Eq. (18)). This segmentation was additionally dilated by a cubic structuring element of size 7 (Eq. (19)) and removed from the approximate white matter segmentation. (Eq. (20)).
Then, a dilation with a cubic structuring element of size 7 is applied. (Eq. (21)). The dilation result is recorded as D 4 . An erosion with a cubic structuring element of size 7 is computed (Eq. (22)). Then, an ultimate volume opening of the result of this erosion step is computed (Eq. (23)). This ultimate erosion step is stored as D 7 . In parallel, we remove D 4 from the segmentation of the whole larvae (Eq. (24)). After this removal, an erosion with a cubic structuring element of size 7 is computed to improve the size of the boundary between in and out regions (Eq. (25)). The result of this last step is summed with two times D 7 Eq.(26). This step creates a map of inside and outside zones. The results of this step is called D 9 . A marker image is created by keeping the local maxima of D 1 and by attributing to these later values stored in D 9 (Eq. (27)). A morphological gradient of D 1 is computed (Eq. (28)). Using this gradient and previously explained markers, a classical watershed is computed (Eq. (29)).

Segmentation of whole larvae
We chose the dorso-ventral orientation, because with lateral acquisitions, the eyes are too opaque. This would result in important signal loss, hiding the data behind them. Automated segmentation of whole larvae was performed as described in the Methods section. These segmentations were compared with the results of the 12 manual segmentations performed by two neuroanatomists. A visual example of these comparisons is shown in Fig. 2 C. Time spent to segment whole larvae was measured for both specialists and for automated segmentation on 6 samples. These times are listed in Table 2. Automated segmentation of these files was almost 3 times faster than performing manual segmentation.
MCC, GB2 and DSC results are shown on Table 3. MCC scores were high enough to consider this algorithm as accurate. However, some regions of the larvae have not been included in the segmented domain, such as the end of the tail (arrow head in Fig. 2.B and C) and, to the contrary, over-segmentation was also observed (arrow in Fig. 2.B and C). However, false positive and negative regions were very small compared to the size of the ground truth segmented domain.  A. Lateral view of a 3D acquisition of a zebrafish 5-dpf larva. B. Automated segmentation of whole larva. Magenta pixels correspond to automatically segmented pixels. This segmentation presents over-segmentation on the head of the fish (arrows) and miss some parts of the tail (arrow head). C. Comparison between a manual segmentation and the automated segmentation. Magenta pixels are automatically segmented pixels, green ones are manually segmented ones. Grey pixels are the overlap between manually and automatically segmented pixels. This comparison confirms the over-segmentation on the head part (arrows) and the missing part on the tail (arrowhead). D. 3D visualisation of the segmentation of whole larva scale bar = 250 µm

Confocal microscopy acquisition and signal decrease in optical sections
In conjunction with the fact that zebrafish larvae are small and relatively transparent, our fast tissue clearing protocol allows acquisition of images with a conventional confocal microscope. Despite this protocol, attenuation of the excitation light as well as the emitted light occurred, leading to signal loss in the deepest regions of the specimens, as shown in a transversal X-Z optical section computed from the image pile ( Fig.1.A). This prompted us to develop algorithms improving contrast in the deepest regions of the specimen. To better esti-mate signal loss, we chose to record samples from opposite directions (left and right sides of the larvae). Two successive acquisitions were indeed possible because fluorescent dye (DiI) was not bleaching out.

Median value of signal loss follows an exponential decrease
The overall signal loss vs. depth into a specimen follows a Lambert law (see Fig. 3). Our contrast correction algorithm started at the depth in the sample where the maximum median value was observed. This maximum value is not located on the first layer of the segmented domain because the skin is not strongly labelled by DiI and thus has low luminosity values.
Using this information, we fitted an affine function to the log of median grey value for each depth above maximum median value using least median of square to evaluate if our estimation followed an exponential decrease.

Contrast correction improves registration
Since the opposite images of the same specimen differ only in the direction of the z-axis, we hypothesized that a rigid registration (rotation, translation) should be sufficient to align two opposite acquisitions. Initially, the rigid registration of datasets without correction led to strongly misaligned elements. This induced a duplication of the retina, tecta (Fig.4.A) and spinal cord (data not shown), making any further analysis impossible.
We show in Fig.4 that computing a segmentation-based contrast correction on both acquisitions of a specimen prior to registration led to a visible improvement of the alignment (Fig.5.B), in particular in the retina and tecta (box in Fig.4A).

Deep signals restoration after segmentation-based contrast correction
In this section, we provide another illustration of image improvements following the application of our algorithm. We applied the parameters of registration, which were previously generated from the alignment of the images acquired from the left and right sides to the initial images, and applied the segmentation-based contrast correction algorithm to both side acquisitions. Fig.5 shows that segmentation-based depth-dependent contrast correction strongly enhances the signal-to-noise ratio in elements like retina (arrowhead), tectum (arrow) or ventral part of the brain (Compare C to A and D to B). However, a blurring effect, which is due to the depth-dependent changes the point spread function (PSF), is not corrected by this method. We however estimate that this will not alter subsequent volumetric analysis of the white matter much, because this blurring effect does not affect significantly the global shape of the positive domain.

Segmentation-based contrast correction allows the improvement of manual segmentation
To measure the efficiency of our contrast correction method, one of our neuroanatomists segmented the white matter of the same specimens before and after contrast correction. Results of the difference between manual segmentations are visible in Fig. 6. Clearly, the white matter of the ventralmost part of the brain was under-segmented when no contrast correction was applied. From the shape of the ventralmost part of the segmented white matter, we hypothesize that it belongs to the hypothalamus, the ventralmost part of the brain and therefore, that the whole brain was segmented following contrast correction.
Before contrast correction, a large part of the ventral part of the brain white matter was missed due to contrast loss over depth. This loss is visible in Fig. 6 A and B.
As with the segmentation of whole larvae (refer to section 4.1), two trained specialists manually segmented the white matter on raw data. The times spent to perform manuals and automated segmentation of 6 samples are visible on Table 4. Automated segmentation of white matter on these files is at least 4 times faster than performing manual segmentations. To measure the miss-targeting effect on raw data, DSC, GB2 and MCC were computed, comparing these segmentations to manual ones performed after contrast correction. These measurements are shown on Table 5.
MCC and DSC results show a strong impact of the contrast correction. For example, MCC shows a median discrepancy of 17.76% for the first neuroanatomist. These errors are impacting volume and position

Accurate white matter segmentation
We compared manual segmentations performed on contrast corrected data to corresponding automated white matter segmentations. Visualisation of automated white matter segmentations and their comparisons against manual results is shown on Fig. 7. Manually segmented regions present a large variety of shapes, such as linear elements in the tail (data not shown) or small dots in the dorsalmost part of the brain (arrows in Fig. 7 A and C) which were not well detected by our algorithm. MCC, GB2 and DSC scores of these comparisons are shown in Table 6. Median MCC, GB2 and DSC scores are respectively, 86,18%, 85,66% and 86,11%. By comparing these results with those given in Table 5, we con- Figure 6: Segmentation-based contrast correction improves manual segmentation accuracy. Orange represents the segmentation before contrast correction, blue is segmentation after contrast correction and white is the overlap between these two categories. A. Comparison between manual segmentations performed before and after segmentation-based contrast correction. Ventral part of the brain, deep in the specimen, was missed without contrast correction. B. Ventral view of the surface rendering of the comparison between manual segmentations performed before and after segmentation-based contrast correction. A large part of the ventralmost part of the brain was missed without contrast correction and was present after the application of the contrast correction procedure. scale bar = 250 µm.

Figure 7:
Comparison between manual segmentation performed on contrast corrected data and the automated segmentation including contrast correction. Magenta represents automatically segmented voxels (including contrast correction). Green represents manually segmented voxels. White represents overlap between automated and manual segmentations. Ultimate volume opening removed the nonconnected elements of the brain (arrows). A. Parasagittal optical sections. B. Ventral views of surface renderings. C. Lateral views of surface renderings. scale bar = 250 µm clude that the automated segmentation of the white matter is more accurate than the manual segmentations performed without contrast correction.

Conclusion
In this paper, we describe a simple and accurate method to automatically segment and then measure the shape of an entire 5dpf zebrafish larva and its white matter. This algorithm opens the route to large scale white matter volume analyses in the framework of High-Content/High Throughput Screening. Measurement of these volumes could for example be used to perform toxicological screens for analysing the impact of a drug on the central nervous system or to study the effect of a mutation affecting overall brain morphology.
In addition to this labelling of the white matter, labelling of the grey matter using an anti-HuC antibody labelling could be additionally performed, enabling the analysis of whole brain volumetry.
Also, in the near future, new markers of brain sub-domains may provide refined analysis of potential abnormalities. Finally, our method will provide the ability to parcellate brain regions, in conjunction with both classical statistical, image registration-based atlas methods, as well as learning-based methods thus opening the route to more refined morphometric analyses. All these avenues for future work are currently under investigation. Table 5: Measure of the MCC and DSC for the segmentations of the brain white matter comparing manual segmentations performed on raw data by two neuroanatomists to manual segmentations performed on contrast corrected data by one of these neuroanatomists.