Spatial mapping and quantification of developmental branching morphogenesis

Branching morphogenesis is a key driver of organ development across taxa (for reviews, see Affolter et al., 2009; Andrew and Ewald, 2010; Costantini and Kopan, 2010; Miura, 2008). During this process, epithelial organ primordia elaborate within the context of a surrounding mesenchyme, initiating a complex genetic crosstalk between these populations that drives tissue differentiation. Many genetic pathways have been proposed to shape branching morphogenesis, on the basis of either organ culture experiments or subjective assessments of the morphology of branching structures. However, quantifying the precise contribution of a given gene (or environmental factor) to branching has been constrained principally because of a lack of computational approaches to assess changes in organ shape and structure in three dimensions. Conventionally, ex vivo organ culture is employed to solve this problem because the flattened tissue generated by this approach is amenable to visualisation using conventional microscopy. However, these thin, avascular organs develop outside of the foetal milieu, significantly disrupting normal spatial tissue organisation. Consequently, they often poorly reflect in vivo organ development (Al-Awqati and Goldberg, 1998; Short et al., 2010). As a result, our understanding of the spatial development and patterning of organs, even during normal development, is often poor. Moreover, these limitations affect our capacity to assess how branching morphogenesis is affected by specific genetic perturbations or by alterations in the foetal environment.

Perhaps the most impressive study of branching morphogenesis undertaken to date has been the investigation carried out in the developing mouse lung by Metzger and colleagues (Metzger et al., 2008). In this study, static images of many stained whole-mount foetal organs collected across developmental time were effectively annotated ‘by hand’ to profile the branching of the pulmonary epithelium. This study revealed a remarkable developmental stereotypy regulated in a distinctive and reproducible manner by different genes and genetic programmes. Despite the significance of this study, the application of this kind of ‘manual’ approach is not suitable for the types of analyses required to routinely assess the impact of genetic or environmental insult. More importantly, it does not permit the objective, quantitative measurement of branch characteristics such as curvature, angle and rotation, which are essential for a full appreciation of development in 3D space. Such ‘higher order’ measurements will be of great importance in allowing us to understand fully how an organ develops within the embryo.

To date, attempts to accurately measure aspects of branching morphogenesis have focussed on subjective hand annotation (as described for the lung above), reconstruction and annotation of serial sections (Blanc et al., 2012; Sims-Lucas et al., 2011) and our own initial attempts at computationally deriving branching datasets. The ‘by hand’ approaches are incredibly time consuming requiring significant user input to annotate branching events and our initial computational attempts have been inflexible in their capacity to capture information from the remarkably diverse 3D forms generated by branching embryonic epithelia. In this article, we describe an entirely new method by which these complex measurements can be made. We demonstrate the application of this approach to the study of branching in a number of organ systems, including the lung and kidney, and we employ it to identify previously unrecognised contributions of a gene to organ formation. This advance will provide researchers with the capability to accurately measure and assess branching datasets generated using recent advances in imaging technology. Moreover, it promises to provide crucial insights into how organs normally develop in 3D space and precisely how this remarkably complex process is disrupted either by genetic mutation or through changes to the foetal environment.

Foetal mouse kidneys (n=4) and lungs (n=4) were stained as previously published (Metzger et al., 2008) (Short et al., 2010) using pan-cytokeratin (Sigma) and E-cadherin (Invitrogen) antibodies, respectively. Embryonic day (E) 9 chicken lungs were dissected in PBS, fixed in 4% paraformaldehyde for 1 hour at 4°C, then stained using the same protocol as for mouse lungs. E14.5 mouse salivary glands were processed as for the kidneys and stained for E-cadherin (cadherin 1 – Mouse Genome Informatics). All secondary antibodies were AlexaFluor 555 (Invitrogen) conjugated.

Samples were imaged in a Bioptonics 3001 OPT scanner (Bioptonics, UK), at resolutions appropriate for the size of the organ. OPT datasets are reconstructed using nRecon (Bruker microCT, Belgium). When not rendered using our software, rendering of volumetric data was performed using Drishti (http://anusf.anu.edu.au/Vizlab/drishti/). Circular phylogram data were rendered using JStree (http://lh3lh3.users.sourceforge.net/jstree.shtml) and linear phylogram clusters with TreeGraph2 (http://treegraph.bioinfweb.info/).

Unless otherwise indicated, all statistical tests comparing two groups where a two-sided P-value is noted are the result of a Student’s t-test with Welch’s correction for potential unequal variances.

The TS program is available to academic users upon request. Datasets may be filtered using a 3D median filter to remove localised signal noise, as either a 3³ or a 5³ voxel neighbourhood. Cropping is performed to enhance processing speed and reduce memory requirements. The software performs a watershed fill in xyz dimensions through all greyscale levels above a baseline to produce a filled greyscale volume. Initially, to determine the centre of the tubular volumes, 3D Gaussian blurring is applied and a vector gradient is calculated using an inverse-square-law field. The calculated vectors follow the 3D Gaussian blur intensity gradient and point towards bright regions (representing the centre of filled tubes). The vectors following the gradient have greater magnitude when the change from dark to bright is stronger. This vector field calculates a weighted average of gradient vectors at all points at which the weight obeys a configurable inverse power law. The result is a vector field influenced by the global structure of the sample and shaped by the local minimum size structures in the sample, a value that can be controlled by the user. The vector field is then used to calculate a convergence field, which marks points as being on the centre line or centre plane of the samples and these act as seed points for the skeletonisation. These are iteratively traversed outwards in a 3³ voxel neighbourhood honouring the direction of surrounding field vectors. The result is a thick fully connected object skeleton. The vector field also provides a cost function for thinning this fully connected object skeleton to generate a final binary structure.

This ‘estimated’ binary voxel skeleton can then be edited within the software in a 3D interactive viewer, in which the tree ‘root’ is also manually set. Then the ‘corrected skeleton’ is traversed from the user-defined root point of the tree, identifying branch points and terminal (tip) end points, and a graphical representation of the tree is constructed from the branch/end points (nodes) and segments (edges). Segment IDs are assigned to a continuous sequence of skeletal voxels located between branches and/or terminal points in the tree. A second editor can then be used to curate the segmentation, allowing the user to merge segments together. This is advantageous for analysing segments in organs, such as the lung, in which domain branching generates a shared central bronchial ‘airway’ with offshoots. Further editors allow the user to define body plan axes in the structure (dorsal/ventral/left/right) by orienting the structure on the screen, choosing a reference voxel in the skeleton and executing the command. These data are then propagated throughout the tree and exported with the tree measurements. Similarly, a labelling function allows initial voxels defining segments of a structure (such as the lobes of the lung) to be labelled and this label is similarly allocated to subsequent segments/nodes derived from that point.

The segmented skeleton is overlaid on the filtered scan data, and a star-convex fill is used to propagate outwards from the central skeleton to the boundaries of the local tubular volume, which is defined by the background level intensity used in the cropping process and is limited by the maximum feature diameter setting. This generates a binary volume surrounding the segmented skeleton and voxels are segmented into zones based on the nearest skeletal voxels, which are pseudo-coloured (Fig. 1B) and quantified to provide the total volume of each segment. Splines are then derived from the skeleton voxels using a recursive least squares fit of a spline of a degree proportional to the number of skeletal voxels in the segment scaled for minimum feature size and limited to ten control points. At the terminal points of the tree (where the skeleton ends centrally in the tip), the spline is extended to reach the volume surface boundary (Fig. 1B). Measurement of splines is more accurate than estimating metrics from a coarse binary voxel skeleton, as splines are smooth and measurements are performed at sub-voxel precision. All length, local angle and curvature data are directly derived from the smooth 3D splines, which describe more accurately the ‘biological’ shapes of the branches of such organs. Determination of diameter also relies on the use of splines, where the spline is traced at 0.1 voxel precision along the length of the spine and the diameter detected by testing the distance from the spline to the surface from a normal vector.

Individual angles from the incident are measured for the start direction of each segment (given by the first derivative of the spline). For bifurcations, this is the angle between the start direction of two splines coincident at a branch point (local bifurcation) or relative to the next branch or terminal point (global angle). Where present, dihedral angles are measured between two successive sets of bifurcations. Planes are fitted to the direction vectors of the two splines coincident at a branch point, and then again at the parent bifurcation branch point (local dihedral angle). The dihedral angle is the angle of rotation between these two planes. The same measurement is repeated but using directions to the next branch of terminal point (global dihedral angle). Curvature is measured as the average per-unit-length acceleration of a particle following the trajectory of a spline.

Tree Surveyor generates a summary spreadsheet that contains segment ID notation, segment predecessor (parent segment) ID, segment name, generation (calculated from root), inverse generation (calculated from tip), length, median diameter, volume, curvature, angle (from incident), local bifurcation, global bifurcation, local dihedral, global dihedral, a tag if the segment is a terminal segment, a tag if a segment is an offshoot segment and, for offshoot segments, a global and local offshoot angle. Propagated directional orientation information and labels are included in this output. These metrics are also exported in GraphML format, which includes attributes for graph edges such as the segment ID name, the node IDs the edge is connected to, length, median diameter and volume. Node attributes provided are the location in 3D space (xyz coordinates provided, including those for the termini and extensions of termini to the surface boundary), local and global bifurcation angle, local and global dihedral angle (coincident at the node) and the source and target edge IDs. Basic tree architecture can also be exported in Nexus format so the skeletonised organ plan can be visualised in a flat dendrogram-style format using external tools (see Fig. 4B; supplementary material Fig. S1A and Fig. S2).

User controls enable optimisation for different organs, in particular the ‘minimum feature diameter’ setting, which controls many of the Gaussian and vector gradient stages. This was set to 9 μm for E15.5 kidneys and 15 μm for E14.5 lungs. The application and scoring of the voxel field is controlled by sliders for ‘features’ and ‘seeding’. For the kidney and lung samples, the features slider was set to –3.000 and seeding was set to 0.800.

We set out to design a computational pipeline in which we could accurately assess the complex morphological measures essential for a complete understanding of branching morphogenesis. We elected to use optical projection tomography (OPT) (Sharpe et al., 2002) to examine branching morphogenesis, principally because it is suited to the visualisation of mouse embryos (Alanentalo et al., 2008; Short et al., 2010; Walls et al., 2008). However, the approach is equally applicable to other imaging modalities that can produce isotropic 3D datasets of labelled branching epithelia. The resulting program (termed Tree Surveyor, hereafter TS) is designed to interrogate branching morphogenesis in 3D. The principal challenge in assessing branching is to refine a complex volume into a skeleton from which measurements can be derived. The TS program employs a 3D vector-based skeletonisation approach (loosely based on 2D optical character recognition algorithms) to overcome this (Fig. 1A). OPT datasets of stained epithelia in developing organs are watershed filled across the range of pixel intensities and then Gaussian blurred in 3D. The gradients established by this blurring define vectors that are analysed to delineate a convergent, weighted centre of the branching volume, which in turn defines the skeletal dataset (Fig. 1A). A ‘root’ for the skeleton is manually selected as the basis for subsequent analysis. To the best of our knowledge, no commercially available software can accomplish the refinement of these structures across different organ types, including our previously published, kidney-specific attempt at such analysis (Short et al., 2010). Our approach can be fine-tuned for different tissues by modulating minimum feature size, the strength of the Gaussian field and the weighting placed on the vectors used for skeletonisation. Altering these settings also allows the program to adapt to local differences in staining intensity or expression.

No automated skeletonisation approach is ever 100% accurate. We therefore developed an editor that visualises the generated skeleton relative to the volumetric data used to derive it, permitting elements to be deleted, joined or added (Fig. 1B; supplementary material Movie 1). This ensures that TS generates quality-controlled datasets with absolute accuracy and fidelity. Manual curation is required for <2% of branch tips and only very occasionally for internal skeletal elements. Individual branches of a tree often develop into different functional units in the adult organ (as in lung airways). To study such ‘segmental metrics’, TS defines branch nodes in the skeleton and projects the interconnecting segmental information back onto the original scan data (supplementary material Movie 2), dissecting it into individual volumes (Fig. 1B; supplementary material Movie 3). To refine spatial measures further, curves of best fit (splines) are calculated for each of the segments, providing a more ‘organic’ representation of the skeleton. Direct ‘point to point’ lines between branch points define ‘global’ branch angles, whereas the splines are used to derive initial ‘local’ branch angles (supplementary material Fig. S1B). Splines also determine further characteristics, including segmental curvature and length. Finally, in terminal branches, the splines are projected to the very surface of the source volume (using their terminal direction), providing accurate tip distance measurements (Fig. 1B; supplementary material Movie 3). An initial survey demonstrated that TS was able to successfully skeletonise a number of developing organs in different organisms, including the mouse lung, kidney and salivary gland and the complex air sacs of the foetal chick lung (Fig. 1C-F).

The end result of these analyses is a set of metrics that define the characteristic features of any given branching epithelia. TS produces a spreadsheet with information for all segments, including their hierarchy, generation number, length, median diameter, volume, curvature, local and global bifurcation angles, local and global dihedral angles and whether or not a given branch is terminal. As each segment is assigned a generation number (relative to the root) and inverse generation number (relative to the tip), metrics can be derived for all segments belonging to a specific generation. TS also produces a GraphML file containing all of the data in the spreadsheet as well as the coordinates of segments, branch and tip nodes in 3D space. Finally, it generates a Newick format file which can be used to plot a 2D ‘phylogeny’, in which the length of a given segment in 3D is directly proportional to the length of the rendered 2D branch (Fig. 4B; supplementary material Fig. S1A and Fig. S2).

To determine whether TS could map branching and uncover previously unrecognised features of development, we first examined the branching ureteric bud epithelia of foetal kidneys (at E15.5, n=4 littermates). Cytokeratin-stained organs were OPT scanned and skeletal (Fig. 2A) and segmental (Fig. 2B) metrics assessed using the TS program. Manual editing of the skeletons was required for <2% of individual elements (principally to refine branch tips), highlighting the accuracy of the program in mapping branch elements. In this developmental snapshot, we found that the ureteric tree is composed almost entirely of bifurcations (supplementary material Fig. S1A) and that branching events tend towards orthogonal, with an average rotation of 66.5° (Fig. 2C). The degree of dihedral rotation in the kidney at this developmental time point was remarkably consistent across branch generations (Fig. 2D). During successive generations, branches become shorter (P=0.0002), thinner (P<0.001) and straighter (P=0.0002) with concomitant decreases in volume (P=0.0045) (Fig. 2E,F). The global bifurcation angle is remarkably uniform until the sixth generation, at which point it significantly increases (P=0.02) (Fig. 2F; supplementary material Fig. S1C,D). At the same time, we found no significant change in local branch angle size between generations (one-way ANOVA, P=0.7). The divergence between these measures correlates with increased curvature of the segments (Fig. 2F).

We next applied our program to the analysis of foetal mouse lungs. Unlike the kidney, the lung is composed of continuous segments with multiple lateral branches located along a proximo-distal axis. To take this structural variation into account, we introduced a ‘segment editor’ into TS that allows bifurcations that might ordinarily be considered a ‘branch’ (with independent offspring) to instead be defined as a continuous segment with ‘offshoots’ (supplementary material Movie 2). Owing to the lobar nature of the organ, we implemented a tool to define body axes at any point within the branching structure (supplementary material Movie 2). Subsequent branches inherit the axis, and their direction away from the axis is calculated to provide a picture of the position of branches relative to the ‘organ body’ plan.

To test these refinements, we stained and analysed foetal mouse lungs (at E14.5, n=4 littermates) (Fig. 3A). Analysis demonstrated the decrease in airway diameter from bronchus to distal tubule (∼80 μm to ∼40 μm) (Fig. 3B) as well as the broad range of terminal airway volumes at this developmental time point (mean of ∼1.99×10⁵ μm³) (Fig. 3C). In addition, we note that the right bronchus (∼141.6 μm) is significantly thicker than the left (∼109.8 μm) (P=0.0032) (supplementary material Fig. S4A), paralleling observations in the adult mouse (Ford et al., 2007). We also undertook a complete analysis of total volumes of the different lobes as well as internal and terminal bifurcation angle and offshoot angles (supplementary material Fig. S4B-H) and generated branching data akin to that defined by Metzger and colleagues (Metzger et al., 2008). Analysis of all branch dihedrals (not including domain branches) showed an average rotation between branches that was remarkably similar to the kidney (61.1° versus 66.5°) with a similar distribution (Fig. 3D).

Using a combination of dihedral angles and the body axis information, we mapped the branching programme (Fig. 4A). Analysis of a portion of the lung [cranial lobe, first anterior branch (R.Cr.AL1)] (Fig. 4B; supplementary material Movie 4) broadly matched the pattern and direction of branching previously proposed (Metzger et al., 2008) (Fig. 4C; supplementary material Fig. S3); however, we did note several differences between individual organs from embryonic littermates. When considering the major branches off the R.Cr.AL1 sublobe, the two major initial anterior branching events appear to be conserved in their orientation (supplementary material Fig. S3); however, the number of subsequent dorsal, ventral, anterior and posterior branches are considerably more variable (supplementary material Fig. S3). Within the clusters formed by the initial bronchiolar branching events, the positions and directions of branches are more varied than previously observed by Metzger et al. in relation to the set body axis (supplementary material Fig. S3).

To determine whether TS can elucidate novel findings about the contribution of a gene to the programme of branching morphogenesis, we examined kidneys from mice heterozygous for a mutation in Tgfb2, the gene encoding Tgfβ2. We have previously reported that Tgfb2 haplo-insufficiency causes a reduction in renal branching, resulting in a decrease in tip number, terminal tip length and decreased terminal branch angle (Short et al., 2010). Analysis using the TS program identified these same defects but also highlighted several other crucial dysmorphologies in the terminal segments, including alterations in global (but not local) branch bifurcation angles (P<0.001), global and local dihedral angles (Fig. 5A) (P<0.001), volume (P<0.001) (Fig. 5B), diameter (P=0.03) (Fig. 5C), curvature (P<0.001) (Fig. 5D) and length (P<0.001) (Fig. 5E). We wondered whether these defects might reflect a simple developmental delay and so compared the mutant datasets with wild-type kidneys with a similar number of terminal branches and generations (E14.5). Although this comparison ‘normalised’ many of the altered measures detailed above (for example, tip volume, Fig. 5B), we still demonstrated an increase in global bifurcation angles (Fig. 5A) (P<0.001), a reduction in tip diameter (Fig. 5C) (P<0.001) and straighter terminal segments when compared with organs of a same ‘notional’ developmental stage (Fig. 5D) (P<0.001). When assessing all branch segments per generation, significant straightening (as shown by a lower curvature metric) is observed from the fifth branching generation out to the tips in the Tgfb2^+/– heterozygote kidneys (maximum P<0.01) (Fig. 5F). On the basis of this analysis, we conclude that Tgfb2 haplo-insufficiency causes a generalised delay in renal development which is coupled to changes in the developmental programme that change the spatial morphology of specific elements of the ureteric tree.

Although our capacity to image branching morphogenesis has improved with advances in microscopy, quantifying the complex 3D epithelial networks that result from this process is an ongoing challenge. We report an analytical platform for undertaking such analysis representing a significant technical advance that will enable a proper assessment of the contribution of different genes and genetic pathways to this process. The TS program itself is remarkably suited to the types of high-throughput analyses demanded by the developmental biology community. Curation of appropriately collected imaging data are typically required for <3% of tips and represents a relatively trivial process in the TS-user interface, wherein the computed skeleton and the actual scan data are co-displayed. We estimate that a competent user can generate a full dataset describing a dozen E15.5 foetal kidneys (each with >500 tips) in a period of 2-3 days. To the best of our knowledge, there exists no other open source or commercial software that can perform this analysis, including our previous (and unrelated) program KAP (Short et al., 2010). We have trialled the capacity of TS to skeletonise and analyse branching structures alongside the ‘filament tracing’ tool in Imaris 7.4, which is designed to analyse neural and filament-like structures. In addition to the inability of this and other programs to effectively skeletonise mid-gestational organs such as the kidney and lung, they also fail to measure dihedral rotation, a crucial measure of organ morphology. Owing to their reduced accuracy in whole-organ skeletonisation, this and other commercial programs require extensive editing of skeletons in order to have them effectively reproduce the branching organ from which they are derived. In light of the laborious nature of hand annotation and the lack of quantitative data this produces, TS will therefore significantly enhance the study of branching morphogenesis.

We have employed the Tree Surveyor program to analyse as many different types of branching morphology as possible, including multiple organs from multiple species. To date, the most difficult structure we have trialled is the E9 chick lung, in which the central lumen provided a particular challenge to skeletal traversal. However, by editing the tree we can set a single length across the bronchus, and connect the detected parabronchial ‘offshoots’ with ease. With this minor exception, we have yet to find a particular branch morphology that overly confounds the skeletonisation mechanism we have devised. This is in stark contrast to our previous report, which employed a completely different ‘erosion’ algorithm for the analysis of the kidney. Although this approach had utility in studying renal branching, it turned out to be particularly stymied by shorter and wider branch morphologies such as those of the domain-branching programme found early in lung development. The new method we report here does not appear to have such limitations. We have favoured the use of OPT data in this study because it generates isotropic datasets (which scale uniformly in all directions); these provide the most accurate metrics describing branching. However, there is no reason why other isotropic computed tomography (CT) or magnetic resonance imagining (MRI) datasets with suitable tissue resolution and contrast could not be analysed using TS. Indeed, carefully collected confocal/light sheet microscope datasets in which z-dimensions are appropriately controlled have proven suitable for analysis. The principal consideration in any of these approaches is the provision of images in which suitable contrast and resolution between the tree and the surrounding tissue has been achieved. To this end, the application of new antibody staining methods for large branching biological samples (Alanentalo et al., 2007) will further increase the utility of this approach.

To demonstrate the utility of this approach, we have employed TS to profile several organs for which development is orchestrated by branching morphogenesis. Our analysis of the kidney reveals that the organ is composed of quasi-orthogonal bifurcations in which a daughter branch rotates, on average, ∼66° relative to its parent. This measure is remarkably consistent and highly reproducible from branch generation to generation, indicating that the TS platform will be useful for detecting modest alterations in the branching programme. In our analysis of the kidney, we find little evidence of trifurcations at this developmental stage, despite previous reports using culture-based approaches that propose a significant role for this type of branching during organ development (Watanabe and Costantini, 2004). At this time point, the kidney has a distinctive morphology in which branches change in shape and thickness from the renal medulla to the cortex, marked by a transition around branch generation 6. The examination of local versus global branch angles is particularly informative in this respect, highlighting very different and more curved branch morphology in the medulla, compared with the very straight branches at the cortex.

The landmark study of lung development by Metzger and colleagues (Metzger et al., 2008) profiled a stereotypic programme of organ development that follows a prescribed pattern of lateral domain branches and either planar or orthogonal bifurcations. We have used TS to examine the morphology of the foetal lung at E14.5 to explore this pattern formation further. It has been proposed that bifurcation events are strictly planar and orthogonal in this developing organ. However, in our analysis we did not observe the binomial distribution among bifurcation events that one might expect if branching was strictly one type or the other. It is important to note that, unlike the quasi-temporal study conducted by Metzger and colleagues, we are considering the actual branch rotations in situ. Consequently, the dihedral angles we have measured might not be as they have originally formed across developmental time. Nonetheless, we noted that the degree of rotation between one bifurcation event and the next in the lung was remarkably similar to that of the kidney. The uniformity of orthogonal rotation in the lung and kidney (despite variations in other measures of branching) suggests that a driving space-filling mechanism shapes the spatial relationship between branches in these organs. Our analysis of the overall 3D morphology of specific parts of the foetal lung agrees with the work of Metzger and colleagues at a general level, but we find much greater variation in patterning at a finer scale, even between embryos from the same litter, on an inbred background. Although it is formally possible that this might reflect differences in the genetic background between our two studies, it strongly suggests that the level of stereotypy in lung development is less rigid than was previously thought.

To examine the capacity of TS to identify previously unappreciated alterations in branch morphology, we examined a Tgfb2 haplo-insufficiency model of defective branching morphogenesis in the kidney. In addition to the already reported reductions in branching (Short et al., 2010), we profile a further role for this gene in shaping the arrangement of different branches in 3D space. Tgfβ2 is expressed in the developing renal tubules and has been proposed to play an active role in aiding the survival and differentiation of the metanephric mesenchyme (Plisov et al., 2001). This role is consistent with many of the defects we observe, given the close inter-relationships between the mesenchyme and the ureteric bud (UB) tips. We have also used the TS application to show that many of the defects we observe at E15.5 are characteristic of a developmental delay, i.e. the morphology of the kidney is similar to a less-developed wild-type organ with the same number of tips. Crucially though, there are significant morphological differences that persist despite this ‘temporal correction’. The mechanism by which the molecule mediates these alterations remains to be elucidated. However, this analysis illustrates how TS can be employed to tease out subtle differences in the structure of branching epithelia and to provide a better understanding of how dysmorphology evolves over developmental time.

In this article, we detail a new analytical platform that enables the assessment of metrics that are crucial to understanding branching morphogenesis, in organs of very different structure. Although OPT data was employed in this study, TS can be used to analyse volumetric datasets generated by other imaging modalities (confocal microscopy, MRI, microCT) providing there is sufficient contrast to enable skeletonisation. Despite exponential increases in gene expression profiling over developmental time, there exist few quantitative approaches for assessing how candidate molecules might actually shape organ formation. Moreover, surprisingly little is known about the normal spatial development of branching epithelia in developing organs. The capacity of TS to perform accurate measurements of branching in 3D will, therefore, be essential to understanding both normal development and how genetic mutation or altered foetal environment translates into complex morphological changes.

We acknowledge the generous assistance of Ajay Limaye in data visualisation.