Three-dimensional microscopy and image analysis methodology for mapping and quantification of nuclear positions in tissues with approximate cylindrical geometry

Organogenesis involves extensive and dynamic changes of tissue shape during development. It is associated with complex morphogenetic events that require enormous tissue plasticity and generate a large variety of transient three-dimensional geometries that are achieved by global tissue responses. Nevertheless, such global responses are driven by tight spatio-temporal regulation of the behaviours of individual cells composing these tissues. Therefore, the development of image analysis tools that allow for extraction of quantitative data concerning individual cell behaviours is central to study tissue morphogenesis. There are many image analysis tools available that permit extraction of cell parameters. Unfortunately, the majority are developed for tissues with relatively simple geometries such as flat epithelia. Problems arise when the tissue of interest assumes a more complex three-dimensional geometry. Here, we use the endothelium of the developing zebrafish dorsal aorta as an example of a tissue with cylindrical geometry and describe the image analysis routines developed to extract quantitative data on individual cells in such tissues, as well as the image acquisition and sample preparation methodology. This article is part of the Theo Murphy meeting issue ‘Mechanics of development’.

JV, 0000-0002-8924-732X Organogenesis involves extensive and dynamic changes of tissue shape during development. It is associated with complex morphogenetic events that require enormous tissue plasticity and generate a large variety of transient three-dimensional geometries that are achieved by global tissue responses. Nevertheless, such global responses are driven by tight spatio-temporal regulation of the behaviours of individual cells composing these tissues. Therefore, the development of image analysis tools that allow for extraction of quantitative data concerning individual cell behaviours is central to study tissue morphogenesis. There are many image analysis tools available that permit extraction of cell parameters. Unfortunately, the majority are developed for tissues with relatively simple geometries such as flat epithelia. Problems arise when the tissue of interest assumes a more complex three-dimensional geometry. Here, we use the endothelium of the developing zebrafish dorsal aorta as an example of a tissue with cylindrical geometry and describe the image analysis routines developed to extract quantitative data on individual cells in such tissues, as well as the image acquisition and sample preparation methodology.
This article is part of the Theo Murphy meeting issue 'Mechanics of development'.

Introduction
Tissue morphogenesis is inherent to embryonic development as tissues have to constantly and continuously adapt their shape during this process. This is true early in development during gastrulation when the first embryonic tissues arise and embryo patterning and morphogenesis have just begun (reviewed in [1,2]), as well as later on when the embryo is patterned and new structures and organs are being formed and reshaped as the embryo grows. Established models of tissue morphogenesis include limb bud growth and neural tube or cardiovascular network formation [3][4][5][6][7][8][9]. Nevertheless, considering the principles of morphogenesis are nearly as diverse as the different organ shapes, almost every model studied requires specific tools for analysis. For example, tooth bud formation or optic cup morphogenesis both involve very specific cell behaviours that have not been described in other developmental contexts [10,11]. The plasticity of multicellular tissues manifested as a global response at the tissue level relies on individual cell behaviours. These behaviours are precisely regulated and coordinated in space and time by a combination of genetic and mechanical cues. This requires intra-and inter-cellular communication achieved through biochemical signals or directly via force transmission through the cytoskeleton. In response to internal or external cues, individual cells can actively or passively change their shape and size, swap location (through cellular rearrangements/intercalation) and be added to (through cell division) or removed from tissues (live cell extrusion/apoptosis) (reviewed in [12 -15]). Developing tools allowing the tracking of these changes in real time and in living animals is thus essential for understanding the principles of morphogenesis at the cellular scale. The diversity of structures and organs that are necessary to build a fully functioning living organism is immense and to meet that demand tissues assume a wide variety of three-dimensional (3D) geometries. This is particularly true for the cardiovascular system where the geometries of the vascular tissue network are highly variable in single embryos [7,8,16 -20]. To further advance our understanding of tissue morphogenesis, we must be able to follow individual cell behaviours in tissues with more complex 3D geometries. This has generated a demand for better 3D image analysis tools, which has also been fuelled by the recent advancements in microscopy techniques in general and in live-imaging of whole tissues and organisms in particular [21]. Such demand has led to the publication, in recent years, of several image analysis tools dedicated to the study of morphogenesis in 3D [22][23][24][25][26][27][28]. For example, optical mapping at the cellular level on the developing zebrafish heart during heart-looping stages in zebrafish has demonstrated that regional cell-function specificity is acquired concomitantly with structural patterning during heart maturation [29]. Three-dimensional heart morphogenesis is particularly challenging to follow due to the heartbeat while developing and several approaches from microscopy to image analysis tools have been specifically developed to assess heart morphogenetic changes in living embryos [20,30 -35]. A recurrent challenge has been to transform complex 3D datasets into two-dimensional (2D) projections in order to simplify the analysis. For example, transformation of 3D Cartesian coordinates projected into spherical coordinates and then into 2D coordinates has been used to study retina morphogenesis [36]. Similarly, methodology developed to map cilia position and orientation on the inside surface of the spherically shaped zebrafish left -right organizer has enabled the creation of mathematical models to predict physical limits of flow sensing inside the organ [37,38]. Other approaches have focused on minimizing the perturbations introduced by long-term live-imaging and have developed methods to track cell migration on fixed tissue during follicle morphogenesis in Drosophila, achieved by tagging extracellular matrix proteins followed by image analysis that involves 'flattening' of the tissue from 3D to 2D thus revealing the trails left by migrating cells [39,40]. Other approaches have focused on reducing data dimensionality in order to deal with the processing of large datasets produced by fast-acquisition imaging methods such as light-sheet microscopy [40,41]. Such advances have enabled the study of tissue morphogenesis on organ or embryo surfaces that have arbitrary shapes as well as the disentanglement of signal from different layers [40]. These efforts allowed the merging of data from various samples uncovering stereotypical cell movements that could not have been detected otherwise [41]. However, data analysis often requires tailored approaches and multidisciplinary skills that are not always available in biology laboratories.
Here, we describe the methodology we have been developing specifically to analyse the position of cell nuclei on the surface of a cylindrically shaped tissue: the endothelial cells (ECs) lining the blood vessels. More precisely, we use the zebrafish dorsal aorta (DA) as a model system to study morphogenesis of vascular tissues and the dynamic behaviours of their composing cells. The DA is the first functioning embryonic vessel and, in the zebrafish, it runs along almost the entire length of the embryo body to circulate the arterial blood pumped from the heart. The venous blood is brought back towards the heart through the caudal vein (CV), which is the main venous vessel in the zebrafish embryo and is located contiguous and ventral to the DA (reviewed in [5,6,9]). In the zebrafish, these two main vessels are formed de novo by a process called vasculogenesis. Angioblasts (endothelial precursors) originating in the lateral plate mesoderm migrate to the embryonic midline where they coalesce to eventually form the DA and CV [42][43][44][45].
Our methodology can be of broader interest because many organs/tissues, such as the kidney ducts, airways and intestine, assume or can be reasonably approximated to a cylindrical geometry. In addition, it only requires software that is commonly available and already used by the research community studying tissue morphogenesis. The details of the image analysis methodology are described below as well as the entire procedure for both sample preparation and image acquisition.

Material and methods (a) Embryo staging and maintenance
Fish maintenance and embryo collection were carried out as described in Westerfield [46]. Embryos were raised in E3 medium when kept inside their chorions or in 0.3Â Danieau's buffer after dechorionation and always incubated at 28.58C. To inhibit pigment formation and improve sample transparency and light penetration, the embryos were incubated in a 0.2 mM solution of 1-phenyl-2-thiourea (PTU) in Danieau's buffer. Embryo staging was done according to chronological and morphological criteria as described in Kimmel et al. [47].
All experiments were performed using zebrafish larvae between 28 and 72 h post-fertilization (h.p.f.), which is in compliance with the European directive 2010/63/EU and IGBMC guidelines validated by the regional committee of ethics (Cometh, Illkirch, France).

(c) Sample preparation
Sample preparation was divided into three major steps. First, to improve reproducibility and facilitate embryo mounting for imaging using an upright confocal microscope (see next section 'Image acquisition'), we designed a plastic mould to cast an agarose pattern where the embryos were mounted (figure 1a,b). The blueprint for the mould was designed in AutoCAD (Autodesk) and the mould 3D printed using resin (ASG-PICO-TRAY, Asiga). In particular, the mould design contains small cavities that accommodate the yolk sac allowing the embryos to rest perfectly on their flank, thus ensuring that the anterior -posterior axis is parallel to the bottom of the dish. The size of the cavities rstb.royalsocietypublishing.org Phil. Trans. R. Soc. B 373: 20170332 was designed to fit the yolk sac of 24 h.p.f. larvae when the yolk is largest and most likely to impair a correct mounting, but also works for larvae up to 72 h.p.f. To cast the agarose pattern, we poured approximately 3.5 ml of liquid 1% agarose gel inside a standard 60 Â 15 mm Petri dish (ref. 353004, Falcon) and placed the mould upside down over the agarose gel. Once the agarose solidified, the mould was gently removed (figure 1a(ii)).
In the second step, we selected embryos fluorescently labelled by both transgenes (see section 2b 'Fish lines') under a stereomicroscope equipped with fluorescence illumination (Leica M205 FA) and transferred them into a separate Petri dish using a fire-polished glass Pasteur pipette. Dechorionation was performed by hand when necessary. In brief, two sharp forceps (Dumont #5, Fine Science Tools) were used to pierce small holes on the chorion, which was then gently peeled off without damaging the embryo. Next, we prepared a 0.2 mM solution of PTU in Danieau's buffer containing 0.1 mg ml 21 of tricaine and a 0.7% low melting point (LMP, Invitrogen) agarose solution in Danieau's buffer that was kept liquid inside Eppendorf tubes placed on a dry heating block set to 428C. Next, the embryos were incubated for 5 min in the tricaine solution to be anesthetized and thus facilitate manipulation during sample orientation.
Finally, for the actual embryo mounting procedure, we first pipetted approximately 0.6 ml of the LMP agarose solution inside the mounting chamber of the agarose pattern (figure 1a(ii)), then quickly pipetted in the embryos with a firepolished glass Pasteur pipette using as little volume as possible. The embryos were then oriented, with the help of forceps, such that they lie on their flank with the yolk sac inserted in the small cavities and their heads all pointing in the same direction (figure 1b). Once the LMP agarose was completely set, everything was covered with the tricaine solution to keep the embryos anesthetized during imaging. Note that the tricaine concentration used is just enough to prevent body movements while still allowing normal heartbeat [50].

(d) Image acquisition
For visualization of the DA and the CV, double-transgenic embryos were imaged using an upright Leica TCS SP8 confocal microscope equipped with a Leica 25Â 0.95 NA water immersion objective (figure 1c,d ). Excitation of GFP was done using a multiphoton laser set to 930 nm wavelength (Chameleon Ultra laser, Coherent, Inc.) while mCherry was excited using a continuous   (e) Image analysis Image analysis was made via two major steps. For the first step, we used the commercially available software Imaris (Bitplane) while for the second we used custom-made Matlab (MathWorks) scripts. In the first step, the cell detection algorithm implemented in Imaris (spots function) was first used to automatically determine the 3D position of each EC nuclei in both the DA and the CV, and manual corrections were made as required. In the second step, the Matlab scripts read the 3D nuclei coordinates previously exported from Imaris, to precisely determine the embryonic axes (A -P, D -V and L -R axis) and re-orient the nuclei of each sample to a common reference system to allow for comparison and averaging of the different specimens acquired . Note that each cross-section, size, orientation (angle of the fitting ellipse), ellipticity and centre of mass change significantly along the vessel due to the irregular shape of the vessel along the anterior -posterior axis. The parametric surface allows mapping of the irregular morphology of the DA to a cylinder and average statistics on the different cross-sections along the same vessels and between different samples. This approach therefore overcomes the ambiguities intrinsic in the complexity of the actual geometry, allowing for a systematic analysis.
In the remainder of this section, we describe the macro steps performed with the Matlab scripts in more depth.
First, the DA and CV nuclei positions (3D coordinates) were imported into Matlab (figure 2a). The nuclei positions of the DA were fitted with a least-square cylinder (using the Matlab routine lscylinder.m, LSGE toolbox) (figure 2b). The F-actin channel was then imported into Matlab and several planes, equally spaced (every 4.4 mm; approx. 50 sections: length ¼ 221.43 mm), were defined along the longitudinal (anterior-posterior) axis of the DA and orthogonal to it. The F-actin signal between each plane was then projected on the planes (cross-sections) and used to fit an ellipse (MathWorks file exchange; fit_ellipse by Ohad Gal, based on the Least-Squares criterion) to each cross-section in order to segment the DA (figures 2b and 5b,c). A similar approach can be done using a membrane or a cytoplasmic marker as long as the vessel walls are labelled. Later, the points of the ellipses were interpolated in space with a 3D spline (spline toolbox, Matlab), which provided the parametric surface describing the vessel wall. The resulting interpolated DA was made of several cross-section ellipses that were equally spaced (every approx. 0.5 mm; length ¼ 221.43 mm) along the longitudinal axis (figure 2c).
To define a reference system common for all the embryos  Finally, a common reference system to all the datasets was generated such that the x, y and z axes match, respectively, the dorsal -ventral, the left -right and the anterior -posterior embryonic axes, and the origin of the reference system corresponds to the centre of the DA (figure 2e). Registering all the samples to a common reference system permits clustering and averaging of the data from different samples while avoiding biases and mistakes that are introduced when orienting the specimens in the sample preparation step before imaging (e.g. LMP agarose embedding). Along the same line, our image analysis routine includes normalization of the data by the average segment length, the distance between two consecutive intersomitic vessels (ISVs). The ISV distance was previously obtained for each sample in Imaris using the 'measurement points' function and was then entered in Matlab. The normalization allows for the comparison of different samples that might have different sizes, particularly when comparing control and experimental groups, such as mutants or drug-treated embryos. After the image registration and normalization steps, different parameters can be extracted from the processed datasets. For example, vessel calibre (cross-section area and perimeter) and shape (eccentricity) can be obtained from the ellipses describing the cross-sections of the aorta.
Afterwards, four equally sized regions (dorsal, ventral, left and right) were defined and the DA nuclei were projected on it, and the number and percentage of nuclei in each region could be extracted (figures 2f and 3b). The subdivision into four regions was done such that the arch length of each region is the same on each cross-section and using the most dorsal point of the cross-section as reference. All nuclei positions could be mapped on an average elliptical cross-section (figure 2f ) or on a circle (figure 3a), where each point is well defined by its angular position. The average cell nuclei distribution around the vessels perimeter (cross-section) could therefore be quantified and plotted as a function of the angular position (figure 3c). This procedure is important for comparison between different developmental stages, experimental conditions (e.g. mutants or drug treatment) or to cluster samples of the same stage that may differ in vessel geometry (e.g. calibre).
The described methodology was validated by repeating the analysis on several datasets artificially generated by arbitrary rotations of one of the acquired z-stacks in the AP direction (figure 4a). The artefacts introduced by the analysis were estimated by measuring the distance between the nuclei from the different datasets in the common reference frame (figure 4b), which should be zero. The maximum mismatch observed on our test was of about 40% of the typical dimension of a nucleus (8 mm), which is small compared to the typical perimeter of the DA (60 mm). These results demonstrate that our routine introduces a small level of error that does not significantly affect the final results because the variability introduced by the method is smaller than the biological variability observed (compare plots in figures 3c and 4c).

(i) Instructions to run the Matlab script
Here, we provide detailed instructions on how to run the Matlab script that refer to the critical steps and complement the guidelines provided inside the code. rstb.royalsocietypublishing.org Phil. Trans. R. Soc. B 373: 20170332 The method described above is implemented inside the script '3DcylindricMap.m', which will import the data needed for the analysis (i); segment, interpolate the DA and orient it on a common reference system (ii); analyse the vessel geometry and position of the nuclei in the tissue of interest (e.g. DA), in four equally sized regions (dorsal, ventral, left and right) and plot the number of nuclei as a function of their angular position (iii).
As a preliminary step, download the functions needed to perform the analysis (keep them in a single folder) and create a folder where you will copy the files necessary to run the image analysis code. These include the Matlab or Excel files containing the 3D coordinates of the nuclei in both the tissue of interest (e.g. the DA) as well as in the reference tissue (e.g. the CV in our case), allowing definition of one of the axes (the dorsal -ventral axis in our case). Make sure that, if Excel format is chosen, the coordinates of the nuclei of the tissue of interest are saved in Sheet1 and the coordinates of the nuclei of the reference tissue are saved in Sheet2. Please be consistent in the choice of the file format: you will have to specify if data are imported from Matlab or Excel files and you will not be able to switch while the code is running. In addition, copy a TIFF.file with the z-stack of the marker describing the vessel wall, which will be used to make the ellipse fitting on the cross-sections (e.g. F-actin). An example of the data to run the code is provided together with the Matlab functions. F-actin). Do not forget the file extension; -In 1.5., specify the average distance, in micrometre, between two ISVs to rescale and compare embryos (e.g. L_ISV ¼ 85 in our example) and the TIFF file properties (voxel size: dx, dy and dz); 2. Save '3DcylindricMap.m' and run the code; 3. A window will pop-up in Matlab displaying the middle slice of your F-actin z-stack where you can draw a region of interest (ROI) around the tissue of interest (e.g. DA) by defining sequentially the delimitation points of the region (left click). Once you have drawn the ROI, double click the mouse left button and the ROI will be applied to all the slices in your z-stack (figure 5a). It is important to be precise when defining the ROI, if made too large signal from adjacent vessels (e.g. ISVs) will be included in the ROI and negatively affect the automated ellipse fitting (next step), thus increasing the amount of manual correction necessary; 4. Another window will pop-up displaying an orthogonal view (cross-section) of the F-actin z-stack and a first ellipse will be automatically fitted to the F-actin signal of the tissue of interest, e.g. DA (figure 5b). In the command window, a question will appear, inquiring if you want to correct, type 'y' or 'n' and press enter. If you are not satisfied (figure 5c), manually draw the ellipse by moving the cursor on top of your crosssection and left clicking to define sequentially the points that will be used to fit the ellipse (note that five points at least are needed). This step will be repeated until all the crosssections along the longitudinal axis of the cylindrical tissue of interest are segmented. The number (Npl0) of reference cross-sections is set to 50 by default, but can be changed before running the code, for example, it can be decreased if the geometry of the tissue is regular (section 2.2, where the number of points used to define the ellipses can also be adjusted, Npe).
The script will now run until the end and several results will appear in the workspace (under the names: 'data', 'data2' and 'data3') that contain, respectively, the results for the nuclei position, the cross-section area, perimeter and eccentricity of the individual ellipses as well as for the mean vessel calibre (cross-section area and perimeter) and shape (eccentricity) (see datasheet key in the electronic supplementary material). Plots of the number of nuclei as a function of their angular position and of the percentage of nuclei as a function of their angular position will also appear in figures 1 and 2, respectively. Finally, you will find, inside the folder created at the beginning of the analysis, three figures plotting area, perimeter and eccentricity of the cross-sections in a folder named 'results', as well as a 'matlab.m' file containing all the parameters used in the analysis. Note that fish are assumed to be mounted with the same orientation as described in the Material and method section and figure 1.

Concluding remarks
The protocol described here provides an effective way to track EC nuclei distribution during DA formation. We are continuously improving our method and our goal is to optimize it to include automated cell contour segmentation and detection of cellular rearrangements. We also plan to develop the analysis of time-lapse datasets together with the possibility to study other features like the localization and number of dividing or extruding/apoptotic cells. In principle, this method could be applied generally to other tissues with tubular architecture, provided that the live-imaging is performed with appropriate temporal and spatial resolution. Comparative studies of the nuclei density changes within vascular tissues or in other tubular organs (such as the kidney, lungs, digestive tract) have the potential to provide new insights into the fundamental rules of tissue morphogenesis.
Data accessibility. A compressed file containing the Matlab scripts and functions needed to run the code is provided in the ESM together with a data example. The blueprint for the mould is provided as stl format. Finally, a pdf file with the datasheet key is also provided.