Zebrafish vascular quantification: a tool for quantification of three-dimensional zebrafish cerebrovascular architecture by automated image analysis

Our understanding of the mechanisms governing vascular development relies heavily on animal models. The zebrafish is an excellent model for three dimensional (3D) in vivo vascular imaging (Huisken et al., 2004; Gut et al., 2017) due to advantages such as optical clarity of the embryo, multiple vascular transgenic lines (Chávez et al., 2016; Chico et al., 2008) and advanced microscopy, including light sheet microscopy (LSFM) that can visualise the entire vasculature. Despite its small size, the zebrafish embryo possesses a complex vascular architecture that expands rapidly over time. Some vascular territories, such as the trunk, display a simple and highly stereotyped pattern, but regions such as the cerebral vessels are far more complex, with greater inter-individual and inter-group variability.

Although acquisition of high-resolution imaging of the vasculature in zebrafish and other preclinical models has advanced considerably, objective quantification remains limited. Most studies rely on subjective visual assessment to identify whether an experimental manipulation affects vascular development. If such effects are detected ‘by eye’, quantification may be performed manually on selected vessels, but this biases towards detection of obvious phenotypes. Furthermore, manual quantification of even a small number of vessels in a single embryo is time-consuming, and does not quantify vascular anatomy elsewhere. Vascular phenotypes are commonly reported using representative images of individual animals, rather than comprehensively reporting inter-individual and inter-group variability.

A robust automated approach that co-registers individuals to display a ‘map’ of regions of vascular similarity and variability, and to quantify 3D vascular anatomy would substantially address these limitations. However, to achieve this, many challenges need to be considered. LSFM datasets are often large, rendering data handling, processing and analysis computationally demanding. The cerebral vessels of a 3 dpf zebrafish measure ∼600 µm laterally and the dorsal cerebral vasculature measures ∼280 µm axially. Vessel diameters typically range from 5 to 20 µm, showing either a double-peak (lumenised) or single-peak (unlumenised) cross-sectional intensity distribution with ∼1 µm thick vessel walls (Kugler et al., 2018, 2019a, 2020a). Thus, to provide sufficient information, images are acquired with ∼0.33×0.33×0.6 µm (x,y,z) voxel size and 1920×1920×460 voxels, resulting in images of about 3.5 GB per embryo. File size increases exponentially with additional fluorophores/channels, time-points, experiments, views, etc.

Analysis of preclinical imaging data often lacks assessment of data properties, such as noise, artefacts and motion. A further challenge is that typically only 3-15% of voxels in vascular imaging datasets are vascular, limiting the information available to drive analysis approaches and making processing potentially prone to bias towards background, rather than vascular information (Haixiang et al., 2017; Rahman and Davis, 2013; Tetteh et al., 2019). These issues have led most quantitative studies to analyse 2D projections, an approach that fails to capture the complexity of the 3D vasculature and is confounded by overlying vessels. 3D quantification is clearly preferable, but is technically and computationally more challenging.

Very few previous studies have attempted to quantitatively characterise the zebrafish cranial vasculature. Tam et al. (2012) and Chen et al. (2012) focussed on characterising small vascular regions, instead of the entirety of the brain vasculature, in confocal microscopy images using commercial software. We previously presented a 3D enhancement and segmentation approach for the cranial vasculature in zebrafish transgenic reporter lines (Kugler et al., 2018, 2019a). Daetwyler et al. (2019) presented a method to segment the zebrafish cerebral vasculature using machine learning (dual ResUNet) but did not include further quantification.

Accurate segmentation is the foundation of subsequent analysis and quantification. We therefore previously developed and validated a zebrafish vascular segmentation workflow and examined its sensitivity, robustness and accuracy (Kugler et al., 2020a,b preprint). We applied this to the transgenic Tg(kdrl:HRAS-mCherry)^s916 (Chi et al., 2008) using Sato enhancement and Otsu thresholding (Kugler et al., 2019a; Sato et al., 1998; Otsu, 1979). Our segmentation approach handles lumenised (double-peak cross-sectional intensity distribution) and unlumenised (single-peak cross-sectional intensity distribution) vessels by conversion of double peaks to single peaks, artificially ‘filling’ all vessels, allowing vascular segmentation independent of lumenisation. We confirmed the robustness of our approach with respect to fluorescence signal, sensitivity to true biological differences and ability to extract vessels at different scales.

Building on this validated segmentation approach, we now present the first open-source 3D image analysis pipeline for zebrafish vasculature, called ‘ZVQ’. We apply this to the cerebral vasculature and show that it reliably provides: (1) registration of multiple embryos, allowing examination of vascular patterning similarity/variability between individuals or groups; and (2) global and/or regional quantification of eight vascular parameters (similarity, volume, surface area, density, branching points, length, radius and complexity).

To demonstrate the scientific utility and throughput of our tools, we applied these to datasets from 13 different experiments and quantified the following: (1) changes in vascular architecture during embryonic development from 2 to 5 dpf; (2) the effect of absent blood flow (Sehnert et al., 2002); (3) pharmacological manipulations of components associated with angiogenesis, vasoconstriction and endothelial cell (EC) architecture [i.e. vascular endothelial growth factor (VEGF) inhibition (Nakamura et al., 2006), Notch inhibition (Geling et al., 2002), F-actin polymerisation inhibition (Morton et al., 2000), myosin II inhibition (Kovács et al., 2004), changed osmotic pressure by glucose increase, and altered membrane rigidity by DMSO treatment; (4) the effect of previously published morpholino antisense oligonucleotides targeting jagged1a (jag1a) (Geudens et al., 2010), jagged1b (jag1b) (Geudens et al., 2010), dll4 (Geudens et al., 2010), notch1b (Quillien et al., 2014) and collagen- and calcium-binding EGF domain-containing protein 1 (ccbe1) (Guen et al., 2014); (5) regional vascular topology of the mid- versus hindbrain; and (6) left-right vascular symmetry.

To allow dissemination and wider application of our method, we implemented it as a workflow in the open-source image analysis software Fiji (Schindelin et al., 2012), which contains a custom graphical user interface (GUI), accompanied with comprehensive workflow documentation. In summary, we describe the first comprehensive 3D quantification approach for zebrafish cerebrovascular quantification that generates a far larger quantity of data than manual analysis in a fraction of the time. It detects quantifiable differences of a magnitude that is difficult or impossible to detect by eye, and provides a range of novel tools to analyse and present these data. These tools are likely to aid discovery of biological insights into vascular development.

We established the ZVQ workflow in the open-source image analysis framework, Fiji (Schindelin et al., 2012) (Fig. 1A). Briefly, ZVQ builds on our published optimised segmentation approach using Sato enhancement and Otsu thresholding (Kugler et al., 2019a; Sato et al., 1998; Otsu, 1979) that established the first validated segmentation approach for the zebrafish vasculature. This was applied to all data presented to extract the vasculature.

After segmentation, individual embryos are registered to a template animal (see details below). A user-defined region of interest (ROI) can be specified for analysis of specific regions. Following segmentation and registration, to examine patterning differences between animals or groups of animals, a single population average map (PAM) is generated from all animals in a specified group (Fig. 1A). This PAM shows areas of reproducible and variable vascular patterning within this group. Comparison of the PAMs between experimental groups allows assessment of regional vascular patterning differences (specific examples are shown below).

Following generation of the PAMs, automated quantification of a range of vascular parameters was performed (Fig. 1B). These are: (1) vascular similarity (S_n) – the volume of vascular overlap quantified by Dice coefficient (a measure of inter-animal variability in vascular anatomy in the specified ROI); (2) vascular volume (V_n) – vascular voxels classified using segmentation as vascular and non-vascular; (3) surface area (A_n); (4) density (ρ_n) – the ratio of vascular volume to total volume of interest; (5) network length (L_n) – the number of centre-line voxels, derived by 3D thinning of the segmented vasculature; (6) branching points (BP_n) – points where vessels split up into daughter branches; (7) average radius (R_n) – vascular thickness given by the distance from a local centre-line (or vessel radius midpoint) to corresponding vessel walls; and (8) complexity (C_n) – number of concentric shell intersections from the centre point using Sholl analysis, where the subscript n represents a measurement in the n^th fish of the group under study (Fig. 1; Fig. S1).

Once established, we applied ZVQ to a range of experimental datasets to examine its performance in visualising and quantifying biological differences between groups of embryos. ZVQ builds on our optimised segmentation approach, using Sato enhancement and Otsu thresholding (Kugler et al., 2019a; Sato et al., 1998; Otsu, 1979), that presents the first fully validated segmentation approach for the zebrafish vasculature, and was used for all data in this article.

We used ZVQ to examine changes in cerebral vascular development between 2 and 5 dpf. We first used both manual rigid landmark-based and automatic rigid registration methods to co-register 2, 3, 4 and 5 dpf embryos (six per group, registering five individuals to a single target embryo) (Fig. 2A-C; Figs S2, S3), to generate population average maps (PAMs) of these groups. Examination of these PAMs allows identification of regions of similarity and variability within the group (Fig. S3C-G; Movies 1-4). Quantification of Dice score showed higher inter-group similarity after automatic registration (Fig. 2D; 2 dpf manual P>0.9999, automatic P=0.2154; 3 dpf manual P>0.9999, automatic P<0.0001; 4 dpf manual P>0.9999, automatic P=0.3962; 5 dpf manual P>0.9999, automatic P>0.9999) and a lower coefficient of variation (CoV) after automatic registration (2 dpf, 38.18%; 3 dpf, 26.21%; 4 dpf, 29.23%; 5 dpf, 17.93%) than manual registration (2 dpf, 84.02%; 3 dpf, 81.20%; 4 dpf, 26.21%; 5 dpf, 18.39%).

Unlike examination of individual animals, examining PAMs of groups from 2 to 5 dpf allows identification of regions of conserved vascular patterning. Fig. 2C shows PAMs from six animals each at 2, 3, 4 and 5 dpf. These demonstrate a high degree of vascular similarity at all time points, particularly of perineural vessels, such as the primordial midbrain channel (PMBC), primordial hindbrain channel (PHBC), middle cerebral vein (MCeV), mesencephalic vein (MsV) and central arteries (CtAs). Over time, the expansion of the hindbrain vasculature between 2 and 5 dpf, with an increase in the number of central arteries, could be observed.

We next examined the ability of ZVQ to detect vascular topology parameters, which were expected to increase in volume, length, branching points and complexity in the ROI of the entire brain over 3 days of embryonic development. We found that vascular volume (Fig. 2E; P=0.0008; one-way ANOVA), network length (Fig. 2H; P<0.0001; Kruskal-Wallis test) and branching points (Fig. 2I; P=0.0082; Kruskal-Wallis test), but not surface volume (Fig. 2F; P=0.4885), density (Fig. 2G; P=0.2041), average radius (Fig. 2J; P>0.9999; Kruskal-Wallis test) or complexity (Fig. 2K; P>0.9999; Kruskal-Wallis test), significantly increased from 2 to 5 dpf. This suggested quantification by ZVQ is able to detect biologically meaningful changes in vascular parameters with statistical significance, even in relatively low group sizes of six embryos.

Together, this shows for the first time that landmark- and object-based rigid inter-sample registration can be used to bring the vasculature of different embryos into one spatial coordinate system, with automatic registration performing better than manual landmark-based registration. This ability to register multiple embryos allows the consistent automated placement of regions of interest and produces a PAM that allows inter-group conservation of vascular patterning to be assessed. Furthermore, ZVQ rapidly provides quantitation of multiple parameters that show statistically significant changes during early embryonic development in a manner many times faster than manual quantification.

Our workflow was implemented to allow quantification either without registration (requiring manual ROI selection for each individual fish) or with registration (using one ROI selection per group). We therefore compared the performance of our quantification using both non-registered and registered datasets. This showed good agreement of data before and after registration (Fig. S4A-F), with a reduction in coefficient of variation (CoV) of volume (before registration, 11.14%; after registration, 5.98%), length (before, 12.03%; after, 3.56%) and branching points (before, 14.78%; after, 6.21%), but not surface (before, 5.96%; after, 6.02%), density (before, 3.80%; after, 7.57%) or diameter (before, 1.89%; after, 2.02%).

This suggests that registration results in some minor alteration of the quantified parameters. Conversely, quantification after registration reduces manual work by allowing one ROI per group to be applied to all samples rather than specification of the ROI for each individual.

As 3D registration of the zebrafish brain vasculature has never been performed before, we examined which similarity measurement might be most appropriate to assess registration outcomes in our data. We quantified Jaccard index, Dice coefficient, total overlap, mutual information, sum of squared differences, mean square error and structural similarity, finding similar outcomes in all metrics (Fig. S4G-M).

After showing ZVQ is able to display inter-group differences in vascular anatomy using PAMs and to detect differences in quantitative vascular parameters between different developmental stages, we next examined whether ZVQ can be used to display and quantify biologically relevant differences in embryos of the same developmental stage subjected to experimental perturbation. We first quantified the effect of absent blood flow induced by tnnt2a morpholino antisense knockdown. Fig. 3A shows PAMs generated by co-registering uninjected control (n=15), control morphant (n=18) and tnnt2a morphant embryos (n=10). Visual analysis shows similar PAMs, achieved by 3D image registration and represented as MIPs, between untreated and control morphants. However, the tnnt2a morphant PAM shows substantial anatomic alteration. Although there is high inter-embryo similarity in midbrain vessel patterning and PMBC position in the controls (both injected and uninjected), this is reduced in tnnt2a morphants (blue arrowhead), along with reduced head size (PMBC-to-PMBC distance; Fig. 3A, cyan dotted line). In contrast, the stereotypic patterning of the basilar artery (BA) (green arrowhead) and PHBC (magenta arrowhead) are conserved in the absence of flow. Automated quantification showed that Dice coefficient (Fig. 3B), vascular volume (Fig. 3C), vascular surface volume (Fig. 3D), branching points (Fig. 3F), network length (Fig. 3G), average vessel radius (Fig. 3H) and vascular complexity were statistically significantly decreased in tnnt2a morphants (Fig. 3I). The only parameter not statistically significantly altered was vascular density (Fig. 3E). This confirms that ZVQ detects biologically relevant effect sizes induced by a previously published (but not quantified) physiological manipulation. We next examined the ability of ZVQ to detect the effect of pharmacological treatments. To test this, we performed chemical inhibition of VEGF signalling, which is known to inhibit vascular growth (Bower et al., 2017; Ellis and Hicklin, 2008; Gerhardt et al., 2003). To induce a less pronounced effect on vascular patterning than in the tnnt2a morphants, we used very brief (2 h) VEGF-receptor (VEGFR) inhibition with 250 nM AV951 at 4 dpf. Visual assessment of the PAMs (six/group) showed no obvious difference in vascular patterning or topology between control and AV951-treated groups (Fig. 4A), suggesting that visual assessment would not detect patterning differences. However, after quantifying vascular parameters with ZVQ, vascular volume, surface volume, branch points, network length and average radius were all statistically significantly reduced (Fig. 4; Table S1). Together, this showed that quantification by ZVQ is able to detect inter-group differences that would be challenging to detect visually.

We next quantified the effect of other pharmacological manipulations known to influence vascular development (Fig. 5A; Table S1). We examined the effect of Notch signalling inhibition, which was shown to result in hyper-vascularisation (Bray, 2006; Watson et al., 2013). Visual assessment of PAMs (n=6) showed no overt phenotype (Fig. S5), but quantification showed an increase of almost all measured vascular parameters, but only a statistically significant increase in vascular volume and surface volume (Fig. S5).

We then examined the impact of inhibiting F-actin polymerisation and myosin II cytoskeletal components, which are pivotal for EC stability and the vascular smooth muscle cell (vSMC) contractile machinery (Sacharidou et al., 2012). Visual assessment of PAMs (n=6) suggested, for both treatments, an overall preservation of vessel patterning but with regional alteration in the intra-neural midbrain region. Vessel diameters in the PAMs were also observably reduced in comparison with controls. Indeed, when quantifying the vascular architecture upon inhibition of F-actin polymerisation (latrunculin B; Fig. S6) and myosin II (blebbistatin; Fig. S7), a statistically significant reduction of almost all parameters was found.

Next, we studied the impact of increased osmotic pressure using glucose, which is thought to result in cellular swelling, potentially impacting angiogenesis (Scallan et al., 2010; Severinghaus, 1995), with a previous study in zebrafish showing that only exposure for over 96 h impacted vascular topology in tectal vessels (Chhabria et al., 2018). Visual assessment of PAMs (n=6) showed no overt phenotype. Quantifications confirmed that 24 h glucose exposure did not alter cerebrovascular topology significantly (Fig. S8; 1.31%).

Last, we examined the impact of decreased membrane rigidity (DMSO) on vascular topology as it results in increased membrane permeability (de Ménorval et al., 2012) and in transient osmotic absorption in vessels (Glass et al., 2006) and increases membrane permeability (de Ménorval et al., 2012), which could potentially change at least vessel radius. Again, no discernible phenotype was observed in the PAMs and quantification showed that 24 h exposure did not statistically significantly change vascular topology (Fig. S9).

We next examined the effect of previously published morpholinos on the vascular architecture, to further assess the performance of ZVQ and provide quantitative insights into the effects of these morpholinos on cerebral angiogenesis.

Jagged1 has previously been shown to be a proangiogenic regulator in mice, but so far no study has examined the distinctive roles of the zebrafish orthologues jagged1a or jagged1b in zebrafish angiogenesis. Visual assessment of PAMs of both jagged1a and jagged1b morphants suggested patterning alterations in intra-neural midbrain vessels. Quantification by ZVQ showed morpholino knockdown of jagged1a decreased network length and branching points (Fig. S10), while knockdown of jagged1b reduced vascular volume, surface volume and network length (Fig. S11). We next examined the impact of dll4 knockdown, which was previously shown to result in hyper-vascularisation in the zebrafish trunk (Leslie et al., 2007) and mouse retina (Lobov et al., 2007). PAMs indicated that dll4 morphants had a smaller head size, with patterning changes in the midbrain, but not hindbrain. Our quantification showed a decrease of all vascular parameters, except radius and complexity (Fig. S12).

Loss of notch1b has previously been shown not to significantly impact vascular topology in the zebrafish trunk (Quillien et al., 2014). Similarly, PAMs displayed no observable difference, and quantification showed no statistically significant difference in the cerebrovascular architecture upon notch1b knockdown (Fig. S13).

Previous work showed that ccbe1 is required for lymphangiogenesis and venous differentiation (Hogan et al., 2009), and loss of ccbe1 impacts the murine coronary artery (Bonet et al., 2018). PAMs suggested alterations of the midbrain vasculature, with quantification of vascular topology after ccbe1 knockdown showing a significant reduction of volume, surface volume, branching points and network length (Fig. S14).

To further study the impact of the examined morpholinos and chemicals on vascular topology, we quantified percentage differences in mean values of all of the measured quantitative parameters (Fig. 5A; Table S1; control morphant to morphant, and untreated to treated; see Movies 1-16), whereas percentage differences were used to allow comparison between groups.

Cluster analysis revealed five main clusters, which can be considered as: (1) inhibitors of anti-angiogenic factors (notch1b morpholino and Notch inhibition), (2) factors inducing cellular changes (osmotic pressure and membrane rigidity changes), (3) factors with severe angiogenic defects (dll4 morpholino, actin polymerisation inhibition, tnnt2a morpholino and jagged1b morpholino), and (4) inhibitors of angiogenic factors (myosin II inhibition, jagged1a morpholino, ccbe1 morpholino and VEGF inhibition) (Fig. 5B). Our data showed that ZVQ provides informative population average maps and extracts meaningful parameters to quantify and describe vascular topology.

Next, we examined regional vascular topology to examine the ability of ZVQ to quantify vascular anatomy in smaller, specific regions of interest. Our initial manual measurements of local vascular growth had shown that growth was comparable between fish (Fig. S2), thus embryos could be compared with each other on a scale.

We quantified vascular parameters in the midbrain and hindbrain regions at 3 dpf (Fig. 6A). As this has never been carried out before, it was unclear whether vessel radius, network length, branching points or volume would be similar in these regions. We found vascular network length (P=0.6168; Fig. 6B), branch points (P=0.0717; Fig. 6D) and vascular volume (P=0.0717; Fig. 6E) to be not significantly different, while the average vessel radius was significantly higher in the hindbrain (P=0.0008; Fig. 6C). To understand the reason for this, comparing vessel radii further revealed that the BA was the main contributing factor for the increased average vessel radius in the hindbrain (Figs S15, S17). We used the same approach to quantify the vasculature of the left and right brain at 3 dpf (Fig. 6F). We found no statistically significant difference in network length (P=0.6830; Fig. 6G), average vessel radius (P=0.5665; Fig. 6H), number of junctions (P=0.8816; Fig. 6I) or vascular volume (P=0.5730; Fig. 6J). To examine left-right vascular symmetry in more detail and over time, we then developed an image analysis workflow to assess topological similarity of the left and right vasculature automatically (Fig. 6K-P). This was achieved by mirroring the vasculature from the right brain hemisphere to the left by horizontal transformation (Movie 17) and extracting vessel overlap by voxel majority decisions (Fig. 6Q,R). This suggested that peripheral/perineural vessels show a strong degree of left-right symmetry from 2 to 5 dpf (Fig. 6S). Quantification of vascular volume (Fig. 6T) and network length (Fig. 6U), showed no statistically significant difference from 2 to 5 dpf, although the Dice coefficient between the left and right vasculature was low (Fig. 6V). Thus, our quantification approach allows for novel regional insights of vascular topology, showing that the average vessel radius is higher in hindbrain vessels than midbrain vessels, but that the vasculature is highly symmetric from 2 to 5 dpf.

We also examined whether regional differences between sample groups could be extracted by image calculation of overlap and difference (Fig. S16), finding that regions of similarity and overlap could indeed be calculated. This allows researchers to identify regions of interest, which could be used to examine these extracted regions in more detail via specific ROI selection.

Last, we examined whether information could be extracted from individual vessels. Theoretically, various approaches for this exist. For example, (1) vessel start and endpoints could be extracted, and cubic ROIs drawn automatically (Fig. S17A); however, this is inaccurate because of vessel proximity and curvature. (2) Individual vessels could be extracted using a brain atlas (Fig. S17B), as previously shown for mice (Todorov et al., 2020); however, zebrafish brain atlases are only available for 6 dpf and healthy controls; an alternative could be to use a vessel-based atlas, which is likely to work at least for major vessels, but not for smaller non-stereotypic vessels. (3) Individual vessels could be annotated manually and annotations transferred to other vessels (Fig. S17C). We examined this last option using simple neurite tracer (SNT) (Longair et al., 2011), finding it, in theory, applicable, but highly laborious (annotation for one fish taking ∼8 h; Fig. S17D-I).Together, our approach allows for regional studies of the cerebral vasculature, but future work is needed to allow extraction and annotation of individual vessels.

As our aim was to provide a useful and easy-to-use image analysis approach for the community without prior computational knowledge, the following steps were undertaken to allow uptake by the community:

• 1. ZVQ was implemented as one Macro in the open-source software Fiji, which is free to use and is cross-platform applicable. Additionally, the workflow was designed in a modular fashion, allowing users to run one or multiple analysis steps at a time.

• 2. A graphical user interface (GUI) was designed to allow users to select the required steps using drop-down boxes without the need to code (Fig. S18).

• 3.The code is provided on Github, accompanied by an in-depth step-by-step workflow documentation (explaining individual steps, input and output); importantly, the code will be maintained by the lead author, but is free to be modified and adapted by other researchers, depending on their data and analysis needs.

• 4. Example data are provided to allow users to examine the workflow first, before applying it to their own data.

• 5. We have included reviewer and user feedback to improve the presented workflow for use.

In this work, we present ZVQ – the first biomedical image analysis workflow that allows 3D zebrafish cerebrovascular quantification, including: (1) inter-sample registration and generation of population average maps; (2) quantification of parameters (vascular volume, surface volume, density, network length, branching points, radii and network complexity); and (3) quantification of regional vascular parameters comparing midbrain with hindbrain, left with right and regional vascular topology, using the freely available image analysis software framework Fiji.

Our pipeline builds on a validated segmentation approach, applicable for all examined ages and data. Should users wish to further improve segmentation outcomes, this could be achieved by raw data quality improvement, such as (1) prevention of pigmentation using chemicals or zebrafish non-pigmented lines, or (2) multi-view image acquisition and fusion.

When examining registration methods, we found manual landmark-based as well as automatic rigid registration were applicable to our data, and suggest that automatic registration is preferable based on the decreased CoV and the significant reduction in user interaction. However, as for all automatic approaches, registration outcomes should be visually and, where possible, quantitatively assessed. Future work is needed to examine whether registration outcomes could be improved using iterative registration or bespoke solutions optimised for sparse vascular datasets.

Registration resulted in low Dice scores for both registration methods. However, unlike in the medical field, no reference values were available to examine our registration outcomes, as we here present zebrafish vascular inter-sample registration for the first time. To address this, we validated our segmentation approach with experimental data from the same sample after displacement and examined other similarity measurement scores, including Jaccard index, total overlap, mutual information, sum of squared differences, mean square error and structural similarity. As all of these showed similar results, we opted to display the Dice coefficient, owing to its implementation in the free open-source software Fiji and computational efficiency. We anticipate that data sparsity of the segmented vasculature is one contributing factor for these low Dice scores, as small shifts between vessels caused either by inter-sample variation or residual registration errors will affect the overlap more drastically than would be the case for more continuous structures such as brain tissue. Future work is needed to examine Dice coefficients on vessel-segment levels, to study voxel-wise comparisons, to explore alternative similarity metrics more suited to a sparse vascular network and to develop registration methods that utilise higher order transformations.

Our registration approach allows examination of regional similarity, as well as facilitating the development of a vascular atlas. Future work could use double transgenics to allow other information, such as expression patterns, to be co-registered to the vascular atlas, as has been performed for a zebrafish brain atlas (Randlett et al., 2015; Marquart et al., 2017).

We showed ZVQ is able to extract data from different ages (2-5 dpf examined). We anticipate that the change from 2-3dpf in vascular density is likely to be due to the higher head curvature at 2 dpf, rather than to a smaller volume (Fig. 2E) or head size (Fig. S3). We suggest that the transient dip observed in vascular surface volume at 4 dpf could potentially be due to: (1) the amount of pruning starting at 3 dpf resulting in decreased surface area, and (2) the fact that myogenic responses and vascular regulation start at 4 dpf and might alter surface transiently. However, as there is no ‘gold-standard’ to compare our data against, future work is needed to examine this further.

Our studies examining vascular anatomy in groups of animals with and without experimental perturbation (physiological, chemical or genetic) showed ZVQ is able to detect statistically significant differences with relatively low group sizes and despite relatively mild patterning defects in many cases. For example, we show the first quantitative assessments of decreased and increased cerebral vascular patterning in response to VEGF and Notch inhibition, respectively. This is in agreement with the literature that VEGF is pro-angiogenic, while Notch is angiogenesis limiting (Lobov et al., 2007; Phng and Gerhardt, 2009; Hogan and Schulte-Merker, 2017). As expected, our data show that inhibition of F-actin polymerisation or myosin II both result in a reduction of vascular parameters, and we anticipate this is mainly caused by a mechanical collapse of the cytoskeleton and the vascular lumen.

The unchanged vascular topology upon glucose treatment is in agreement with previous work, which showed 24 h exposure has no significant effect on vascular topology (Chhabria et al., 2018). Additionally, the unchanged vascular topology upon DMSO treatment suggests that DMSO indeed only leads to transient vascular changes, as previously suggested (Glass et al., 2006).

These findings clearly demonstrate that quantification of 3D vascular parameters enables meaningful insights into the effects of a range of manipulations on vascular anatomy, and is sensitive enough to extract true biological differences, resulting in a decrease or increase of the vasculature. When examining morpholino data, we show for the first time that loss of jagged1a and jagged1b results in mild cerebrovascular defects. As jagged1a and jagged1b are paralogues found in zebrafish due to partial genome duplication (Meyer and Schartl, 1999), we anticipate that their simultaneous loss might result in a stronger effect.

Finding that dll4 knockdown results in a reduction of the cerebral vasculature is in contrast to previous work in the zebrafish trunk (Leslie et al., 2007). This suggests that dll4 might have cerebrovascular-specific roles that have been previously overlooked and require future investigation.

Our finding that loss of notch1b does not significantly impact vascular topology suggests that notch1b plays a conserved role in the zebrafish head and trunk vasculature (Quillien et al., 2014).

Having found that loss of ccbe1 results in a decrease of vascular volume, surface volume, branching points and network length. We postulate that ccbe1 has a conserved function in different vascular beds, as previous work has shown that ccbe1 is required for zebrafish lymphangiogenesis and venous differentiation (Hogan et al., 2009), as well as murine coronary artery development (Bonet et al., 2018).

Moreover, we developed an automated approach to study the vasculature in defined regions, illustrating the usefulness of the approach in the midbrain and hindbrain, as well as for comparison of the left and right vasculature. However, the implementation of our analysis approach is not limited to this, allowing the user to select a variety of regions of interest.

Future studies are needed to examine vascular symmetry between the left and right vasculature, as well as to examine stochasticity rules of vascular patterning in more depth, to assess local differences and how these might relate to structural and/or functional brain asymmetry, which is established at 4 dpf (Dreosti et al., 2014). As the zebrafish brain vasculature is an enclosed system, which is in contrast to the human cerebrovascular tree, localised effects, such as by a stroke, are likely to translate to global outcomes. We anticipate that our quantification approach will be feasible to allow examination of left-right symmetry in older zebrafish, which may be more likely to show vascular asymmetry. Our quantification was performed at levels of the entire brain or in subregions, and was sufficient to detect significant differences during development and in response to manipulations. Future work might examine vascular parameters in older zebrafish, other transgenic lines and vessel-specific quantification, which will require an in-depth annotation.

We presented a proof-of-principle example on how to examine individual vessels, but future work is needed to establish brain and vascular atlases for zebrafish at varying ages and of varying vascular phenotypes to allow automatic annotation and extraction on vessels. Currently, atlases for the zebrafish brain exist for only healthy animals at 6 dpf (Randlett et al., 2015; Kunst et al., 2019) and adults (Kenney et al., 2021; Ullmann et al., 2010). This means any analysis in other ages or treatment groups would require creation of new atlases to allow for brain/morphological differences with dedicated experiments and analysis. Alternatively, new atlases could possibly be developed using available information and modelling age-specific changes. However, this would require extensive a priori information and validation. To allow application of the analysis approach by other researchers, we have implemented a GUI, and we provide code, example data and step-by-step workflow documentation.

We describe ZVQ – the first comprehensive workflow for assessing zebrafish 3D cerebrovascular similarity, topology and regional variability. Our data show that robust image analysis is needed to allow objective and quantitative insights. Our results emphasise that the cerebral vasculature is highly complex and that many mechanisms regulate its formation and patterning.

Experiments performed at the University of Sheffield conformed to UK Home Office regulations and were performed under Home Office Project Licence 70/8588 held by T.J.A.C. Maintenance of adult zebrafish in the fish facilities was conducted according to previously described husbandry standard protocols at 28°C, with a 14:10 h (h) light:dark cycle (Westerfield, 1993; Aleström et al., 2020). Embryos, obtained from controlled pair- or group-mating, were incubated in E3 buffer (5 mM NaCl, 0.17 mM KCl, 0.33 mM CaCl₂ and 0.33 mM MgSO₄) with or without Methylene Blue. All experiments were conducted with the transgenic reporter line Tg(kdrl:HRAS-mCherry)^s916 (Chi et al., 2008).

Data described in this section were obtained from Kugler et al. (2019b). Briefly, the impact of nine chemical components and six morpholinos (Genetools) on the vascular architecture was examined. Morpholino injections were performed at the one-cell stage using Phenol Red as an injection marker.

The conditions studied, the compounds used and the n values were as follows. VEGF inhibition was achieved using 250 nM AV951 for 2 h (Nakamura et al., 2006) at 4 dpf (control, n=22; treated, n=23; three experimental repeats). Notch inhibition was achieved using 50 µm DAPT for 12 h (Geling et al., 2002) at 4 dpf (control, n=24; treated, n=24; three experimental repeats). Myosin II inhibition was achieved using 25 µM blebbistatin for 1 h (Kovács et al., 2004) at 4 dpf (control, n=21; treated, n=23; 3 experimental repeats). F-actin polymerisation inhibition was achieved with 100 nM Latrunculin B for 1 h (Morton et al., 2000) at 4 dpf (control, n=13; treated, n=12; two experimental repeats). Osmotic pressure was studied using 40 nM glucose for 24 h at 4 dpf (control, n=18; treated, n=18; two experimental repeats). Membrane rigidity was assessed using 0.25% DMSO for 24 h at 4 dpf (control, n=24; treated, n=22; three experimental repeats).

The following morpholinos were used: cardiac troponin T2a (tnnt2a; 1.58 ng; Genetools; 5′-CATGTTTGCTCTGATCTGACACGCA-3′) (Sehnert et al., 2002); jagged1a (jag1a; 0.1 ng; Genetools; 5′-GTCTGTCTGTGTGTCTGTCGCTGTG-3′) (Geudens et al., 2010); jagged1b (jag1b; 0.8 ng; Genetools; 5′-CTGAACTCCGTCGCAGAATCATGCC-3′) (Geudens et al., 2010); delta-like ligand 4 (dll4; 3 ng; Genetools; 5′-GAGAAAGGTGAGCCAAGCTGCCATG-3′) (Geudens et al., 2010); notch1b (0.25 ng; Genetools; 5′-GTTCCTCCGGTTACCTGGCATACAG-3′) (Quillien et al., 2014); and calcium-binding EGF-like domain 1 (ccbe1) (5 ng; Genetools; 5′-CGGGTAGATCATTTCAGACACTCTG-3′) (Guen et al., 2014). Control MO injection was performed according to the same protocol, with final concentrations as above (Genetools; 5′-CCTCTTACCTCAGTTATTTATA-3′).

Experimental repeats and n values for each morpholino were as follows. Cardiac troponin T2a: 3 dpf uninjected control, n=17; control MO, n=18; MO, n=15; two experimental repeats. jag1a: 3 dpf uninjected control, n=11; control MO, n=12; MO, n=10; two experimental repeats. jag1b: 3 dpf uninjected control, n=21; control MO, n=21; MO, n=21; three experimental repeats. dll4: 3 dpf uninjected control, n=23; control MO, n=23; MO, n=23; three experimental repeats. notch1b: 3 dpf uninjected control, n=21; control MO, n=21; MO, n=21; three experimental repeats. ccbe1: 3 dpf uninjected control, n=22; control MO, n=22; MO, n=22; three experimental repeats.

Datasets were obtained using a Zeiss Z.1 light sheet microscope with a Plan-Apochromat 20×/1.0 Corr nd=1.38 objective and a scientific complementary metal-oxide semiconductor (sCMOS) detection unit. Data were acquired with activated pivot scan, dual-sided illumination and online fusion. Properties of acquired data are as follows: 0.7× zoom, 16 bit image depth, 1920×1920 pixel (∼0.33×0.33 µm) image size and minimum z-stack interval (∼0.5 µm), 561 nm laser, LP560 and LP585. Zebrafish embryos were embedded in 2% LMP-agarose containing 200 mg/l tricaine (MS-222, Sigma). The image acquisition chamber was filled with E3 plus tricaine (200 mg/ml) and maintained at 28°C.

Images were analysed using open-source software Fiji (Schindelin et al., 2012). Data analysis was performed without blinding.

Data were pre-processed by performing a .czi to .tiff conversion and by generation of maximum intensity projections (MIPs) using the Bio-Formats plug-in. Intra-stack motion correction was performed based on scale invariant feature transformation (SIFT) (Lowe, 1999, 2004) using the Fiji plug-in ‘Linear Stack Alignment with SIFT’ with parameters as described previously (Kugler et al., 2018, 2019a).

Vascular enhancement was performed assuming local vessel tubularity using the Fiji Sato enhancement filter ‘Analyse>Tubeness’ (Sato et al., 1998) with an optimised parameter set as described previously (Kugler et al., 2018, 2019a,b). Enhanced images were segmented using Otsu thresholding (Otsu, 1979) as described in previously (Kugler et al., 2018, 2019a).

The following measurements were performed manually using the Fiji line region of interest (ROI) tool and ‘Analyse>Measure’. Growth expansion of the primordial midbrain channel (PMBC) was measured at the inner vascular edge at the position behind the eyes. Distance of the middle metencephalic central artery (MMCtA) was measured at the most anterior point. The MMCtA angle was measured in the right MMCtA from its posterior to its anterior end. Posterior metencephalic central artery (PMCtA) angle was measured from its origin in the centre of the embryo to its right lateral end.

Visual comparison of samples was used to identify and extract 11 anatomical landmarks that were found in all samples and distributed along the anterior-posterior, left-right and dorso-ventral body axis. These 11 anatomical landmarks (Fig. 2) were then used for manual landmark-based rigid registration using the Fiji plug-in ‘Name Landmarks and Register’. The automatic rigid inter-sample registration method, by Johannes Schindelin and Mark Longair, based on the virtual insect brain protocol (VIB) (Schmid, 2010; Jenett et al., 2006), was used with the following parameters: no initial transformation, five best matching orientations for further optimisation, image downsampling two to six times, no ROI reduction (i.e. bounding box selection) for optimisation, using Euclidean measure of difference and showing transformed image.

Target embryos, used as a template for registration, were selected based on (1) sample body-axis orientation along the image axis (anterior-posterior along image y-axis, coronal plane along image z-axis and x-axis); (2) all major vessels visualised in image; and (3) no obvious abnormalities (i.e. vascular patterning visually normal, no haemorrhaging and no gross defects).

To allow comparison between sample groups, all embryos were registered to a template from the control group, while ROIs for quantification were drawn for one sample/template per group. To validate registration, we constructed a new dataset by repetitive image acquisition of the same sample after manual displacement. This allowed quantification of registration outcomes following registration to the initial acquisition. To further examine 3D topological similarity between samples, six embryos were registered to a target and averages extracted.

µ is the local mean and σ is the standard deviation.

The ‘Analyse Skeleton’ plug-in in Fiji (analyse>skeleton>analyse skeleton 2D/3D; implemented by Ignacio Arganda-Carreras (Arganda-Carreras et al., 2010) was used to identify and measure the number of branching points (analyse>skeleton>summarise skeleton). To quantify vessel radii (R_n), EDMs were multiplied with extracted skeletons, resulting in a 1D representation of vessel radii as represented by intensity of voxels.

Intra-sample symmetry was assessed using the segmented vasculature to allow topological rather than intensity-based assessments. After manual selection of the central body axis [centrally along the basal aorta (BA), between L-R posterior communicating segments (PCS), between L-R MMCtA and the bifurcation point of L-R anterior cerebral vein (ACeV)], automatic image rotation was performed to align the anterior-posterior axis of samples with the image y-axis to allow for mirroring. Afterwards, the right vascular volume was automatically mirrored by horizontal image transformation and similarity was derived via voxel overlap.

Semi-automatic tracing and annotation of cerebral vessels was performed using the Fiji plug-in Simple Neurite Tracer (SNT) (Longair et al., 2011). To allow detection of these vessels in other fish (following segmentation, registration and skeletonisation), traces were extracted as vessel-specific ROIs (multiple ROIs to describe x, y and z position). These vessel-specific ROIs were then used to produce rectangular ROIs (x,y; 30×130 voxels) to allow the detection of centrelines in other fish. The ROIs were then used to detect and quantify vessel intensities (from EDM and skeleton) to quantify vessel-specific diameters.

The code is available on Github (https://github.com/ElisabethKugler/ZFVascularQuantification) and is accompanied by in-depth step-by-step workflow documentation (SDoc1; questions and bug reports to be sent to E.C.K.). A doi has been assigned with Zenodo (https://doi.org/10.5281/zenodo.3978278). The workflow is designed in a modular fashion, allowing users to run individual or multiple steps at once, depending on their data and analysis needs. It was implemented in an open-source framework to allow adaption and changes by the community. Example data are available on zenodo (10.5281/zenodo.4108660; https://zenodo.org/record/4108660#.X47XU9BKizc).

To visualise data, maximum intensity projections (MIPs) were generated and intensity inversion was applied as appropriate to give the clearest rendering of relevant structures.

Normality of data was tested using D'Agostino-Pearson omnibus test. Statistical analysis of normally distributed data was performed using one-way ANOVA to compare multiple groups or Student's t-test to compare two groups. Non-normally distributed data were analysed with a Kruskal–Wallis test to compare multiple groups, or a Mann–Whitney test to compare two groups. Analysis was performed in GraphPad Prism Version 7. P values are indicated as follows: *P<0.05, **P<0.01, ***P<0.001 and ****P<0.0001. Based on the mean and s.d. of the data in the examined groups in the assays used, post-hoc power calculations have shown that these assays have at least 80% power to detect a 30% effect difference between groups when group sizes are 12/group (α=0.05). Samples were excluded from analysis if image quality did not allow for reliable quantification. Animals were allocated to treatment groups randomly without selection. Imaging and data analysis were performed unblinded to treatment allocation, often because the effect of treatment was easily deduced from the appearance of the micrograph. To overcome subjective bias, objective automated image analysis was applied where possible. Data are mean±s.d., if not otherwise stated. Cluster analysis was performed using ClustVis (Metsalu and Vilo, 2015). Settings used were as follows: clustering distance – correlation, clustering method – average, tree ordering – tightest cluster first. Clustering validation was performed using relative clustering validation (varying the number of clusters from one to three; one displayed). Image representation was performed using Inkscape (https://www.inkscape.org).

We are grateful to Steven Renshaw and Freek van Eeden for their helpful feedback on the manuscript. We thank Robert N. Wilkinson, Stefan Schulte-Merker, Emily Noël, the Perak lab in Sheffield and Henry Roehl for sharing chemical compounds. The authors thank the Bateson Centre aquarium facility for zebrafish husbandry. We thank the reviewers for their in-depth feedback to improve this manuscript, and the biomedical image analysis as well as ImageJ/Fiji community for valuable feedback. Zeiss Z1 light-sheet microscope was funded via British Heart Foundation Infrastructure Award awarded to T.J.A.C.

The peer review history is available online at https://journals.biologists.com/dev/article-lookup/doi/10.1242/dev.199720.