Development of an ImageJ-based method for analysing the developing zebrafish vasculature

Zebrafish with fluorescently labelled blood vessels provide an excellent model for studying angiogenesis. Most commonly the growth of the intersegmental blood vessels is investigated in response to compounds or manipulation of gene expression and analysed using manual methods, typically scoring the connectivity of these blood vessels to the dorsal longitudinal anastomotic vessel.  Such methods are laborious and best suited to time points after the connectivity of these vessels have been established.  By contrast, reported image processing-based methods are difficult to implement and often depend on specialist software. This study aimed to develop and evaluate a computational method using the freely available ImageJ software to quantify the development of intersegmental blood vessel formation. This methodology developed allowed rapid analysis of vascular development.  The outputs of total vessel length and number of junctions best documented defective vascular development at differing levels of severity and gave comparable results to the frequently used manual approach of calculating percentage connectivity of intersegmental vessels. This ImageJ-based method allowed objective quantitation of vascular network formation in zebrafish enabling a free, straightforward and rapid approach to determine the effect of novel compounds or genetic manipulation of the angiogenic process.


Introduction
The zebrafish ( Danio rerio) is a powerful and widely used model organism for the study of vascular development and angiogenesis. Features such as the embryo's rapid and external development, small size, optical transparency, fluorescent tagging of vessels and methods to modify its gene expression all contribute to its utility in studying vessel development [ 1]. This vertebrate model shares both anatomical features and similar molecular mechanisms of vessel development observed in higher vertebrates making this a valuable tool for understanding gene function in both normal human development and angiogenesis dependent pathologies [ 2,3].
The widespread use of zebrafish in studying the effects of genes or compounds in vessel development Simms Victoria A. Bicknell Roy Heath Victoria L. @ + @ is accompanied by a variety of methods used to analyse images of the zebrafish vasculature. Typically researchers focus on the growth of intersegmental vessels (ISVs) and the formation of the dorsal longitudinal anastomotic vessel (DLAV) in zebrafish, as the ISVs represent the first angiogenic vessels formed in this organism [ 1]. The majority of ISV analysis methods involve manually scoring images of zebrafish embryos, typically counting the number of normal ISVs formed [4][5][6] or determining the percentage of embryos with defective ISVs [ 7]. These methods are laborious and reflect the difficulties associated with analysing complex images of whole organisms. A number of groups have reported software-based approaches to analyse vascular development to facilitate hit identification in high-throughput screening [8][9][10][11]; these include the use of algorithms designed to analyse neurite outgrowth [ 10] and artificial intelligence-based image analyses [ 11]. Software based approaches have also been applied to high resolution zebrafish images to determine filopodia lengths and sprout volume [ 12]. While these allow quantitation of vessel parameters, they require specialist computer programming expertise and/or are each associated with financial cost in terms of the software required.
The goal of this study was to develop a freely available ImageJ-based method to analyse angiogenic vessel network formation in developing zebrafish embryos. Transgenic TG [ fli1:EGFP] fish embryos which have vascular-specific expression of enhanced green fluorescent protein (EGFP) throughout development [ 13] were used in conjunction with confocal fluorescence microscopy to generate high resolution vascular images, which were then subjected to processing and analysis using the Analyze Skeleton 2D/3D plugin of ImageJ. This computational analysis method was tested by using images generated via a developmental time course as well as images of embryos with varying levels of vascular disruption produced either pharmacologically or by translation blocking morpholino oligonucleotides. Of the generated analysis outputs, total vessel length and junction number best allowed differentiation between the varying phenotypes. This method allowed objective quantitation of vascular network formation in zebrafish enabling a more rapid and systematic approach to determine the effect of novel compounds or genetic manipulation of the angiogenic process.
Morpholino oligonucleotides (Gene Tools, LLC [USA]) were reconstituted into distilled water and were microinjected into the yolk of the zebrafish embryos at the 1 to 4 cell stage of development at concentrations of 0.2-0.8 ng per embryo, with phenol red co-injected as a tracer. Embryos were treated with 0.003% 1-phenyl-2-thiourea (PTU) at 24 hours post fertilization (hpf) to prevent pigmentation.

Morpholino oligonucleotide sequences
All morpholino oligonucleotides (MOs) were purchased from Gene Tools. The RNA antisense morpholino oligonucleotide sequences were designed to block the binding of the translation initiation complex and were designed by Gene Tools, unless otherwise noted.

Imaging
Zebrafish embryos were imaged in lateral orientation on the Zen 780 Zeiss confocal microscope loaded onto MatTek dishes. The 10x objective was used with optical sections taken at 3 μm intervals. Optical sections of the zebrafish embryos were compressed using the Zen Zeiss software to generate maximum projection images prior to analysis.
The Fiji (Fiji is Just ImageJ) version of ImageJ which features the plugins Skeletonize (2D/3D) and Analyze Skeleton (2D/3D) [ 16] was used to analyse the vessel network. The polygon selection tool was used to outline vessels of interest within the image, a line was drawn along the dorsal aorta (DA) and around the image to create a region which included the DLAV and other vessels branching from the ISV as indicated (Fig 1 ). An ImageJ macro was recorded which automatically processed and analysed the selected intersegmental blood vessels (ISVs) and the dorsal longitudinal anastomotic vessel (DLAV) within the image. This macro used the following tools in sequence: Clear Outside, Make Binary, Fill holes, Despeckle, Skeletonize (2D/3D) and finally Analyze Skeleton (2D/3D), checking 'show detailed info'. The script for the macro is as follows: run("Clear Outside"); run("Make Binary"); run("Fill Holes"); run("Despeckle"); run("Skeletonize"); run("Analyze Skeleton (2D/3D)", "prune=none show");

Figure 1: Steps of the ImageJ analysis method applied to vascular images of fli1-GFP zebrafish
The sequence of image processing is shown Figure 1 . In the output two tables are generated: a results table (an example is given in Figure 2A ) giving detailed information for each segment and a table of branch information (an example in Figure 2B ). This Analyze Skeleton (2D/3D) tool can be applied to 3-dimensional imaging and it gives outputs in voxels, the 3-dimensional equivalent of a pixel. Since the images analysed in this study are 2-dimensional images we have noted lengths of vessels and branched in pixels in our data figures. A segment is a set of connected vessels, and for each segment the Analyze Skeleton (2D/3D) function calculated the number of the following: branches (slab segments connecting end-points, end-points and junctions or junctions and junctions), junctions (with junction voxels merged), end voxels (<2 neighbours), junction voxels (>2 neighbours), slab voxels (=2 neighbours), triple and quadruple points (junctions with exactly 3 and 4 branches, respectively), average branch length and maximum branch length. Example images showing the analysis nomenclature are shown in Figure 2C . The results were collated in Microsoft Excel and the following parameters were recorded for each image. The segment number was given by the number of rows in the results table; the number of junctions and end points were given by the sum of the junction and end voxel columns, respectively. The total vessel length was derived from the branch information table which gave the co-ordinates and length of every branch; it was calculated by summing the branch lengths. The mean vessel lengths were calculated by dividing the total vessel length by the total number of branches, given by the sum of values in the number of branches column in the results table. 10 embryos were imaged and analysed at each time-point in the developmental time-course study. 5-15 zebrafish embryos were imaged and analysed for each set of morpholino injected zebrafish to knock down vascular related genes.
To compare the ImageJ analysis method with the manual analysis approach, the same zebrafish images from the morpholino injected zebrafish images were analysed by manually counting the number of fully formed ISVs which made connections to the DLAV. The percentage of connected ISVs for each embryo was calculated. Figure 2: Outputs of the ImageJ analysis method to quantify the zebrafish vasculature.
A.An example table of analysis outputs containing the information on vessel segments. B.An example table of branch information C.The region highlighted in the red box in Figure 1 was enlarged to show the identification of the analysis output nomenclature of slab voxels (orange), end point voxels (blue), junction voxels (magenta) and vessel branches (white arrows). Three types of vessel branch are indicated, they occur between end voxels, junction voxels or an end voxel and a junction voxel.

Statistical analyses
All data results were expressed as mean ± SEM represented as error bars. Statistical comparisons were performed using GraphPad Prism 6.0. To determine if there were significant differences between the morpholino and control zebrafish vasculature image analysis results, a Mann-Whitney test was performed. Pearson's correlation coefficients were calculated to determine the significance between a number of different parameters as described below in the Results. Differences were considered to be statistically significant when the p value was less than 0.05. Statistical significance has been indicated by the following: *** for a p value of less than 0.001, ** for a p value between 0.001 and 0.01 and * for a p value between 0.01 and 0.05.

Results
In order to test the functionality and applicability of this ImageJ analysis method to the different stages of ISV development, it was first tested on a series of images which tracked vascular development from initial sprouting of the ISVs at 24 hours post fertilisation (hpf) to the complete establishment of the DLAV at 72 hpf ( Fig 3A ). The images were analysed and the number of vascular segments, end points, junctions, total vessel length and mean vessel length determined. For each time point the mean values of these parameters were plotted (Fig 3B ). From 24 to 30 hpf there was a rapid increase in all parameters except mean vessel length which reflected the rapid overall growth and development of the organism. Over this period a rapid increase in total vessel length is accompanied by a rise in the number of branches and this results in a more modest increase in mean vessel length (the total length divided by the branch number). There was a decline in segment number and end points between 36 to 48 hpf consistent with the greater connectivity of the vasculature. A segment as calculated by this method is defined as a set of connected vessels; this is distinct from an intersegmental vessel, which is a vessel that grows between the somites and collectively these are considered to be the first angiogenic vessels to form in this organism [ 1]. The total vessel length and junction number increased rapidly initially with a slower increase observed from 36 -72 hpf. Finally, the mean vessel length rose more steadily between 24 and 36 hpf, with little change thereafter. The graphs generated for these various parameters using this ImageJ analysis method reflect the visual changes observed in the vasculature, and so the method was next tested on images of fish with varying levels of vascular disruption caused by morpholino-induced knockdown of a number of genes known to play a role in vessel development. A.Fli1-GFP transgenic zebrafish embryos were imaged over a 48 hour time period using confocal microscopy. B.The vasculature was analysed using the ImageJ based approach to quantify the number of vessel segments, junctions and end points as well as the mean and total vessel lengths over time and the mean values plotted (n=10, error bars represent SEM) Scale bar: 500 µm.
In order to generate zebrafish embryos with varying levels of impaired vascular development, a number of different genes with known involvement in this process as described in Table 1were targeted using translation blocking morpholinos. The morpholinos disrupted vessel formation to varying extents, ranging from a more severe phenotype caused by the VEGFA and ROBO4 targeted morpholinos, to a moderate phenotype induced by blocking translation of ELTD1 and ETSRP and minimal effect of blocking CLEC14A at 48 hpf. Representative images are shown in Figure 4A of the morpholino injected embryos alongside embryos injected with control morpholinos. Images from both the control and gene-targeted morpholinos were analysed using the ImageJ analysis method and manual scoring for the percentage of connected ISVs. The manual scoring method showed statistically significant differences for knockdown of VEGFA, ROBO4, ELTD1 and ETSRP but not CLEC14A ( Fig 4B ). According to this scoring, ROBO4 morpholino knockdown gave the lowest percentage of connected ISVs, followed by VEGFA and ELTD1 knockdown, then ETSRP knockdown and CLEC14A knockdown disrupting connection of ISVs the least. The different parameters calculated by this ImageJ analysis method are graphically represented in Figure 4B . The end point and segment number parameters gave only statistically significant differences between the control and targeted morpholinos where the phenotype was more severe: only VEGFA for end points and both VEGFA and ROBO4 for segment number. The total and mean vessel length measurements mirrored the data for the percentage of connected ISVs, with significant differences between control and targeted morpholinos for all genes targeted except CLEC14A. The number of junctions was the most sensitive measure of disrupted vascular development with the greatest differences between control and targeted morpholinos, and statistically significant differences for all the genes examined, including a small but significant difference in the vasculature of the CLEC14A morpholino treated embryos. From these analyses the total vessel length and junction number were chosen as the key parameters to use in further analyses undertaken, both because of their sensitivity and because data from these parameters mirrored the percentage ISV connectivity.

Table 1
A summary of the genes targeted during this study.

ROBO4
The predominant roundabout receptor expressed in the vasculature. Found within angioblasts, in the DA, PCV and ISVs in zebrafish, ROBO4 targeted morpholino caused misdirected and truncated ISVs [ 26] VEGF-A Key pro-angiogenic factor, caused severe disruption of ISVs when knocked down in zebrafish [ 27] Figure 4: Analyses of embryos treated with CLEC14A, ETSRP, ELTD1, ROBO4 and VEGFA specific and control morpholinos comparing the ImageJ analysis method and manual analysis of ISV connectivity.
Morpholino injection at the one-to-four cell stage of embryonic development in fli1-GFP zebrafish was performed to silence CLEC14A, ETSRP, ELTD1, ROBO4 and VEGFA with translation blocking morpholinos or their respective control morpholinos. A.Embryos were imaged at 48 hours post fertilization and representative images are shown. Scale: 500 μm. B.Images were analysed using both the ImageJ analysis method and by manually scoring and enumerating the % ISV connectivity (indicated by the hatched shading). The data obtained is represented in bar charts. Significant differences were determined using a Mann-Whitney test, n=5-15, and error bars represent the SEM.
Morpholino knockdown of CLEC14A and ECSCR, genes known to function in vascular development, resulted in evident defects in vascular development at 30 hpf, but not at the later time-point of 48 hpf when the previous analyses were performed (Fig 4A and 5A ). Analysing the 30 hpf images with the ImageJ method revealed statistically significant differences in both the total vessel length and junction number between control and translation blocking morpholinos for both genes (Fig 5B ). The percentage of connected ISVs was also enumerated, but this analysis was more problematic at earlier time points due to the variability in connectivity of ISVs in the control embryos. A significant difference was observed for CLEC14A, but not ECSCR in this parameter. These data suggest that the ImageJ based approach may have wider applicability to the analysis of the vasculature of embryos at different time points in development, in particular before the ISVs have connected with the DLAV.

Figure 5: Analyses of embryos treated with CLEC14A and ECSCR specific morpholinos at 30 hpf comparing the ImageJ analysis method and manual analysis of ISV connectivity.
Morpholino injection at the one-to-four cell stage of embryonic development in fli1-GFP zebrafish was performed to silence CLEC14A and ECSCR with translation blocking morpholinos or their respective control morpholinos. A.Embryos were imaged at 30 hpf and representative images are shown. Scale: 500 μm. B.Images were analysed using both the ImageJ analysis method to give total vessel length and junction number as well as by manually scoring and enumerating the % ISV connectivity. The data obtained is represented in bar charts. Significant differences were determined using a Mann-Whitney test, n=5-15, and error bars represent the SEM.
Finally the ImageJ analysis method was applied to images of zebrafish embryos treated with sunitinib, the VEGFR2 kinase inhibitor. Embryos were treated with sunitinib at 0.1 and 1.0μM and imaged at 48 hpf. A dose-dependent response was observed with minor defects observed with 0.1μM sunitinib treatment and severely disrupted vascular development was observed with 1.0 μM sunitinib (Fig 6A ).
This level of disruption was reflected by the measurements of total vessel length and junction numbers, with slight reductions in both parameters at the lower dose, statistically significant only for total vessel length and much larger statistically significant reductions in both parameters at the higher dose of 1.0μM (Fig 6B ).  In order to further compare different outputs from the ImageJ analysis method, and to compare these with the manual method, the correlation between the different parameters was assessed. For each of the different studies described above including the developmental time course, the morpholino knockdown and sunitinib treatment; the total vessel length, junction number and segment number were each compared (Fig 7 ). Data points were taken from each time point of the developmental series, the morpholino knockdown treatment (controls were not included) and the sunitinib treated embryos and the degree of correlation between all of these pairs of parameters was determined to be very high, with the highest found between total vessel length and junction number (R=0.984; Fig 7 ).
For the morpholino knockdown data, the Pearson correlation was calculated between the manual analysis giving the percentage connectivity and total vessel length ( Fig 8A ) or junction number ( Fig  8B ) were calculated as R=0.863 and R=0.916, respectively, with statistical significance observed for the correlation between junction number and % connectivity of ISVs.

Figure 7: The strongest correlation was identified between the parameters of total vessel length and the number of junctions.
Correlation graphs comparing the parameters of total vessel length, junction number and segment number were plotted using the combined mean data sets from the ImageJ based analysis using the time course developmental data points (blue), the morpholino knockdown experimental data points (red) and sunitinib treatment (green). The Pearson Product Moment Correlation Coefficient (R) from a two tailed t-test were calculated and statistical differences have been shown to be highly significant, for all pairings shown p ≤ 0.0001. A.A correlation graph showing the total vessel length and the number of junctions from the combined data sets, R=0.984. B.A correlation graph showing the number of segments and the number of junctions from the combined data sets, R= 0.931. C.A correlation graph showing the total vessel length and the number of segments from the combined data sets, R= 0.915. Figure 8: Correlation of the manual analysis calculating the percentage of connected intersegmental blood with the computational ImageJ analysis method parameters of total vessel length and the number of junctions.

Discussion and Conclusions
The aim of this study was to develop and evaluate a computer based method that could be applied to quantify and analyse vascular development in zebrafish with fluorescently marked vessels; in particular we aimed to evaluate the formation of the ISVs as they are the first angiogenic vessels to form in the zebrafish [ 1]. The use of ImageJ software was selected as it is available free of charge to the scientific community and thus available to all.
In this study confocal imaging of zebrafish at various stages of development was undertaken as well as various methods to disrupt vascular development. The images generated were subjected to ImageJ based processing focussing on ISV sprouting and connectivity. Once the vessels of interest were selected the image was binarised, and then the fill holes and despeckle functions were applied prior to skeletonisation. The fill holes function was utilised in this method to remove the additional connections formed within the DLAV after 36 hpf reducing the complexity of this part of the vasculature and better enabling the analysis method to target the ISVs and their connection to the DLAV. The despeckle function was applied to reduce levels of noise. The Analyze Skeleton function generated a detailed output and initially the number of segments, number of end points, number of junctions, total vessel length and the mean vessel length were selected to evaluate as potentially useful parameters to reflect the developing fish vasculature. The analysis procedure was first applied to images of embryos taken over a time course covering the initial growth of the ISVs and their fusion to form the DLAV. Very rapid changes in the developing vasculature occur between 24 and 30 hpf and this was reflected in the rapid increase in value of end points, segment number, junction number and total vessel length over this time period. At 30 hpf many of the connections with the DLAV had been made, and in agreement with other studies we found that primary intersegmental vessel formation was completed at approximately 36 hpf [ 17]. After this time the vasculature continued to develop with vessels such as the intercostal and parachordal vessels forming from 36 hpf. This continued, but slower, development was reflected by a more gradual increase in the total vessel length, junction number and mean vessel length between 30 and 72 hpf. There was a decrease in both the number of segments and end points observed between 30 and 48 hpf associated with the increasing levels of connectivity between the developing vessels during this time.
The use of gene specific morpholino oligonucleotides to block either protein translation or RNA splicing has made the zebrafish a powerful tool with which to probe the function of genes in vessel formation in vivo. We sought to test this ImageJ analysis method of analysing the developing zebrafish vasculature on images from zebrafish treated with translation blocking morpholinos targeting a number of different transcripts which gave rise to varying levels of disrupted vessel formation. Inhibiting VEGFA and ROBO4 caused the most disrupted ISV sprouting, while blocking ELTD1 and ETSRP gave a milder phenotype with inhibition of CLEC14A resulting in little discernible effect at 48 hpf. Images were analysed both by using the ImageJ analysis method and by manually scoring ISV connectivity. There was no difference in the number of end points for the knockdown of any gene except for VEGFA which caused severe disruption, this parameter is thus not a sensitive indicator and useful only to document severe absence of ISVs or potentially to enumerate extensive hyper-sprouting of the vasculature. The differential in segment number between the control and translation blocking morpholinos was evident for all genes except CLEC14A, but significant differences registered only for the more disrupted phenotypes, making this an insensitive parameter with which to evaluate milder defects. By contrast, total vessel length, mean vessel length and junction number were significantly different in the more mildly affected embryos compared to their respective controls. In addition, the differential between the control and translation blocking morpholinos enabled clear ranking of the severity. Junction number was the most reduced parameter when comparing control and targeted morpholino knockdown, with a small but significant difference observed in the CLEC14A treated embryos. The parameters of junction number and total vessel length were enumerated from images of embryos treated with different doses of the VEGFR2 inhibitor sunitinib, and though defects caused by the lower dose of 0.1 μM sunitinib were mild, a significant difference in the total vessel length was observed. A clear relationship between vessel growth and branching was demonstrated by a very strong, almost linear, correlation between these numerically unrelated parameters observed between the vasculature of embryos at different developmental stages and after treatment with the different morpholinos.
In order to compare this ImageJ based method with commonly used analyses of zebrafish vessel development we scored ISV connectivity, a method used in many studies [4][5][6][7]18]. Thus each ISV was examined to determine if it had reached and connected with the DLAV, and the percentage of connected vessels calculated for each embryo. Using this method, significant differences were observed for knockdown of all genes except CLEC14A at 48 hpf, and in this regard results were most comparable to the total vessel length and mean vessel length. In terms of assessing the severity of the phenotype, the overall ranking of the five knockdowns was similar. However, the ImageJ analysis method showed a greater defect in terms of every parameter for VEGFA compared with ROBO4 silencing, while the scoring of the ISV connectivity showed the converse. This highlights a difficulty associated with this scoring system such that it does not reflect the extent of ISV sprouting, failing to differentiate between a lack of sprouting (observed with many ISVs with VEGFA knockdown) or developed ISVs unable connect with the DLAV (as occurred with ROBO4 knockdown). Such difficulties can be overcome with more sophisticated scoring systems (eg as used in [ 19]) but these remain time consuming and laborious. An additional difficulty we observed with scoring ISV connectivity was in respect to the timing of this analysis, at 48 hpf the vast majority of control ISVs had connected with the DLAV while those with induced vascular defects had not, thus making this a useful time window for comparison. We observed that for two genes, ECSCR and CLEC14A, silencing caused transient defects which by 48 hpf were not evident, but significantly reduced total vessel length and junction number were recorded for both silenced genes at 30 hpf. At this earlier time point, ISV connectivity of the controls was more variable and a significant difference in the percentage of connected ISVs was recorded only for CLEC14A knockdown, therefore the ImageJ analysis method is useful for vessel analysis at different stages of early embryonic development.
Use of zebrafish as a tool for screening genetic mutations and toxicants on vascular development has prompted the development of high throughput image analysis approaches for quantifying the vasculature [ 8,9,11,20]. Though these reported analysis frameworks often give excellent visual representations of the analysis they can be in practice hard to implement or adjust due to the necessity to use programming languages such as C++ [ 8]. Some methods use microscopy software to assist with manually measuring the lengths of the vessels [ 20] and such methods which focus on ISV vessel lengths did not generate data on vessel junctions and hence give no indication of vessel connectivity [ 9]. The strategy described in this study is not fully automated and requires user input to select the region or vessels to be analysed; however the image processing was then relatively quick with each image analysed in less than 1 minute and generating an output which can be manipulated to give the parameters of interest.
In summary, an ImageJ based analysis method was developed to quantify vascular development in zebrafish and was compared with manual scoring of ISV connectivity. This enabled a rapid analysis of the developing vasculature with total vessel length and junction number providing the best parameters with which to determine disruption in vessel formation and in particular in ISV sprouting, and giving comparable results to the traditionally used manual scoring method.