3D analysis of the whole subcutaneous adipose tissue reveals a complex spatial network of interconnected lobules with heterogeneous browning ability

Adipose tissue, as the main energy storage organ and through its endocrine activity, is interconnected with all physiological functions. It plays a fundamental role in energy homeostasis and in the development of metabolic disorders. Up to now, this tissue has been analysed as a pool of different cell types with very little attention paid to the organization and putative partitioning of cells. Considering the absence of a complete picture of the intimate architecture of this large soft tissue, we developed a method that combines tissue clearing, acquisition of autofluorescence or lectin signals by confocal microscopy, segmentation procedures based on contrast enhancement, and a new semi-automatic image analysis process, allowing accurate and quantitative characterization of the whole 3D fat pad organization. This approach revealed the unexpected anatomic complexity of the murine subcutaneous fat pad. Although the classical picture of adipose tissue corresponds to a superposition of simple and small ellipsoidal lobules of adipose cells separated by mesenchymal spans, our results show that segmented lobules display complex 3D poly-lobular shapes. Despite differences in shape and size, the number of these poly-lobular subunits is similar from one fat pad to another. Finally, investigation of the relationships of these subunits between each other revealed a never-described organization in two clusters with distinct molecular signatures and specific vascular and sympathetic nerve densities correlating with different browning abilities. This innovative procedure reveals that subcutaneous adipose tissue exhibits a subtle functional heterogeneity with partitioned areas, and opens new perspectives towards understanding its functioning and plasticity.


Details of the segmentation procedure
The stitched image mosaic of an entire tissue acquisition can be as large as 12000x12000x300 pixels, leading to images of several tens of gigabytes. Each elementary acquisition volume ranges between 256x256x300 and 1024x1024x300 voxels and many hundreds are assembled within a mosaic. For this reason, the automatic treatment of the image is much more easily handled with a customized and scriptable home-made C++ code, which we developed using the Cimg library (http://cimg.sourceforge.net/index.shtml).
We first segmented the two tissue regions where there was no "lobular organization" as opposed to the one where lobules could be segmented 1 . A mean 3D filter from Fiji software was set up, whose size could vary depending on the sample. The mean 3D filter was chosen to be significantly larger in the x and y directions than in the z direction to take account of the anisotropy of the tissue in the z direction. Then, a simple thresholding of the graylevel was performed. Since the signal was much lower in peripheral regions than in the core region it was accentuated by using the 3D filter. After this thresholdingsegmentation step, the smallest connected components were removed. Finally a convexhull algorithm using the 3D mesh voxelizer program binvox (http://www.patrickmin.com/binvox/) was used to delineate the core region inside which the lobule segmentation was to be performed.
In order to objectively identify the organization of the 3D lobules, we designed a specific work-flow consisting of the eight steps outlined in Figure 1 a) The stitching was performed with free software 2,3 with a 10% overlap between the adjacent volume of the mosaic and based upon maximum correlation. A pyramid of images with three hierarchical levels was also created so as to save memory allocation in some parts of the work-flow.
b) This step was decomposed into b1 and b2. In b1) the filtering step was adapted.
Because of the large heterogeneity of the image texture within each lobule (see Figure 1a) which, in many cases, is as large as the amplitude of the background noise, there was a need to find an adapted filter which, on one the hand, would smooth graylevels inside the lobules, while, on the other hand, preserving or even improving the contrast at their edges, so as to help their further segmentation. For these reasons, step b1) is crucial, with the additional requirement that image inhomogeneity associated with the background illumination is dealt with. Various methods such as variational or level-set methods and texture filtering could have been used in this first step 4 . We used an adapted Kuwahara texture filter, which allowed local edges to be preserved in the smoothing of lobules 5,6 . For each pixel, this approach considers a current cubic window of size L= L x x L y x L z , centred on the current voxel. This window is composed of N x x N y x N z sub-windows, of size S x x S y x S z (in voxels). On each sub-window, the average grayscale texture is computed and then the central voxel is set to the value of the median of all means. For algorithmic optimality, when performing the computation in one pixel, the benefit of the previous computations is kept for its neighbour. Thus only the contribution of each new surface within each cubic-window is evaluated, whilst disregarding the contribution of lost surfaces, as is usual in sliding window methods. The algorithmic cost is thus O(N x N y N z L 3 ) with N x, N y and N z corresponding to the image dimensions in the x, y, and z directions. Varying the size of the window, the number of sub-windows and their sizes showed that N x x N y x N z = 3x3x1, S x x S y x S z = 3x3x3, so L = L x x L y x L z = 9x9x3 gave the best results for edge preservations and graylevel smoothing within the lobules. In b2) the image was thresholded for further use in the work-flow. Various threshold values were tested. c) A simple threshold was used in order to binarize the filtered lobules ( Figure 1c). d) Then the Euclidean distance map to the binarized lobules was computed so as to be able to separate the various entities by a simple thresholding on the distance map ( Figure 1d1). This procedure is similar to a mathematical morphology closure but permits elongated throats and bridges to be cut while preserving the shape of the lobules.
Step ⁄ for two lobules having labels i and j). The graph-merging procedure consisted of merging two edges when their edge weight was larger than a threshold. This threshold was selected as a quartile of all edge weights. The result is illustrated Figure 1i where the small isolated red watershed region inside the blue region visible in Figure 1h has been merged. Similar ideas of merging areas by a graph representation have already been used in various contexts 7,8 . The final result of the lobule segmentation seems to be improved by this last step, the main interest of which is to eliminate over-connected segmented lobules. The gradient procedure that defines the walls for the watershed might still have some "holes" separating entities that would not be disconnected otherwise.
In the following Let us finally discuss the limitation of the study associated with the lobule segmentation procedure. It must be pointed out that, since several arbitrary parameters (such as the binary thresholds) were used through the work-flow, the results might be different if other parameters were chosen. Due to the lack of benchmark and/or ground truth for segmented lobules, image analysis validation used an expert's visual segmentation. Another limitation of the method concerns the previously described border artefacts of the watershed, which are generic. Since the segmented structures are complex and heterogeneous, the watershed can over-connect some of them, which will thus not stop at the precise frontier of the "real" lobule. The segmentation of the "core-region", the filtering, and the graph merging were implemented to reduce the spread of the watershed, but some artefacts might still persist, resulting in spurious, hopefully limited, very small subunits.

Robustness of segmentation method:
In this section we investigated the influence of parameter variations in the segmentation results. Focusing on fat pad n°1, we modified the arbitrary parameters at steps C, E and I (used in table S1) for three values (5, 10 and 20% variations) and studied their effects at a given slice in the tissue sample. Once a parameter of a given step was changed, all other parameters were kept identical to study its specific influence. Size of the Kuwahara-like mask and its sub-windows as well as the type quantity computed in each sub-window at step B1 were not investigated on this section because of its computational time on the full 3D image and the many parameters involved. However, its qualitative effects were studied on small crops of tissue for several sub-windows sizes, sub-window numbers and quantity computed on each sub-window (variance, mean and median). The previously given parameters were seen to provide the best homogenization of gray levels as well as the preservation of contrasts between entities. The size influence of the smallest connected component was not investigated since it is in direct relation with the percentage of merged lobules in step I. Indeed, keeping larger connected components in step F (respectively smaller) will lead to a smaller (respectively larger), percentage of merged lobules in step I.
Finally, threshold in B2 wasn't investigated as it is mostly a visualization parameter and its influence would be mainly focussed on merging step I. Considering the image dynamic, the relative differences in terms of normalized contact-surface between pre-segmented entities wouldn't change much from the threshold at step B2. Figure S1 shows the segmentation obtained for fat pad n°1 with different parameters at depth z = 487.05 µm. Reference (obtained with parameters described in table S1 and illustrated in figure S2A,  The results show that the most critical arbitrary parameter is the threshold at step B. With a variation of 5%, the final segmentation keeps the main characteristics of the reference segmentation with a large central pink subunit (Figure S1 b1) but a higher threshold by a 20% variation will lead to over segmented subunits as illustrated in Figure S1 b2. Size of gap filling is quite robust and leads to almost identical segmentations ( Figure S1 c1 and   c2) as is the threshold on step E ( Figure S1 d1 and d2). Most of the segmented subunits are the same as in reference segmentation although one main subunit close to the lymph node in reference segmentation (dark blue color) is now divided in two subunits ( Figure S1 d1 and d2). Finally, changing the number of merged subunits ( Figure S1 e1 and e2) lead to having two main central subunits instead of the main pink one.
We also quantify for each new 2D crop image, the relative difference with the reference from computing the relative mutual information (relative-MI) 9 as the ratio of the maximum MI, i.e. the MI between the reference image and itself, and the MI between the new segmented image. The results are displayed on supplementary table S3. The higher the relative-MI is, the closer is the image to the reference. As illustrated in Figure S1, the lowest relative-MI value is obtained for 20% threshold variations at step C (Relative-MI = 82%, Supplementary table S3: Consequences of changes in the value of segmentation parameters. Last line reports the relative mutual information (MI) between the reference segmentation image and the new segmentation image obtained with a given set of parameters (Fig. S5). The relative-MI was expressed as the ratio between the maximum-MI obtained with the two images.

Timing of the segmentation method:
Supplementary

Additional discussion about the State of the art
Image segmentation is a long standing topic in image analysis, of particular relevance to biomedical imaging communities. The workflow proposed in this article is original as it has been dedicated to the specific task at end but follows conventional pipelines. Deeplearning methods disrupt the field of image segmentation because they permit robust, performant, and fully automated workflows when possible ( 10-13 ). Deep-learning methods outperform traditional approaches when it is possible to train the algorithms on huge databases onto which ground truth segmentation is available. Nevertheless these supervised approaches are often a distressing step on large 3D images. Manual segmentation for building training data-set remains a big challenge for specific and dedicated task. When the segmentation concerns cancer tumors detection, it is easy to find among the radiologists community the required man power and motivation to build such ground-truth data-set, because this task is indeed a crucial daily-based, intensely sensitive issue for a huge medical community. Regarding the topic of interest here, there is no chance of building the required training data-base of 3D fat-pad lobules. The building of a training data base is not only a costly and intensive issue, but also requires a huge amount of 3D images not possible to provide in our case. Hence for the issue at end in this very specific set of images, we did not found possible to use a deep-learning approach for 3D-lobules segmentation.