Branching morphogenesis in the developing kidney is governed by rules that pattern the ureteric tree

The increase in organ size associated with the evolution of complex body plans presents challenges for the organism in mediating the efficient exchange of nutrients and the removal of waste. To effectively ‘plumb’ tissues, animals have evolved a developmental mechanism known as branching morphogenesis that serves to break up large tissues to maintain their effective surface area to volume ratio. The branching process typically initiates through the formation of bud-like organ anlage, which invades (and ramifies within) an associated mesenchymal field to form complex networks of epithelial tubes that ultimately allow for cellular exchange. Branching morphogenesis is pervasive in higher eukaryotes and is crucial for the development of organs such as the kidney, lung and mammary gland. In the metanephric kidney, organogenesis begins with the formation of the ureteric bud (UB), which appears as an outgrowth of the nephric duct at E10.5 in mice and during the fifth week of gestation in humans (Saxen, 1987). The UB then grows to contact and invade a surrounding field of cells known as the metanephric mesenchyme, which in later development condenses around the tips of the ureteric tree to form discrete cellular fields known as the cap mesenchyme (CM). The UB tips require the presence of CM cells to undergo ongoing branching (Self et al., 2006). In turn, the CM responds to tip signals to control the self-renewal and commitment of precursor populations to form the nephrons (Kreidberg et al., 1993), which will ultimately integrate with the ureteric tree to form the urine collecting system of the adult organ. Reciprocal epithelial-mesenchymal exchange of this nature is a feature of almost all developmental branching morphogenesis and is central to the functional differentiation of cells in both compartments.

A central question in the study of branching morphogenesis in the kidney is the degree to which stereotypy characterizes elaboration of the epithelium. Do epithelial branching events occur at predetermined places and at defined developmental stages to form a morphologically distinct structure? A landmark study of the developing mouse lung describes a highly predictable program that features three morphologically distinct branch types that are reproducibly employed in a program which establishes a precise pulmonary airway architecture (Metzger et al., 2008). Although this suggests a hard-wired genetic control of the branching process, other reports have argued for a more adaptive mechanism emerging from pulmonary epithelial-mesenchyme interactions (Blanc et al., 2012). Mathematical models have been developed based on such known interactions to explain the process of tip bifurcation in the lung (Menshykau et al., 2012) and kidney (Menshykau and Iber, 2013), and the induction of the ureteric bud (Lawson and Flegg, 2016). These models, which are based on a Schnakenburg-like Turing patterning arising from competition for mesenchymal signal, have so far not been able to explain the broad patterning of the epithelial trees. In the kidney, analysis of epithelial ureteric branching in culture suggests self-avoidance in the growing tips – observations that support an adaptive branching mechanism operating in the organ (Davies et al., 2014). This is consistent with our previous studies that have shown that although branching of the ureteric tree is not spatially stereotypical, its elaboration (based on hierarchy) is indeed highly reproducible (Short et al., 2014). That is to say, while specific branches do not form in the same 3D space at the same developmental time, the extent of ramification of the ureteric tree is highly similar between one organ and the next. These observations suggest that a constrained program of patterning drives the elaboration of the UB.

In order to understand how mouse ureteric tree patterning is repeatable, we have combined three-dimensional imaging, spatial quantification and mathematical modelling to develop a class of state-based mathematical models that describe the bifurcation and elaboration of ureteric trees. We then compared our models that best describe normal branching against ureteric trees from mice carrying mutations in genes thought to be crucial to shaping the branching process. Although highly stereotypical, the self-similar nature of the branching pattern we describe indicates that branching occurs under the influence of local rules, which is consistent with an emergent pattern arising from interaction between neighbouring epithelial and mesenchyme niches. We also show that the characteristic asymmetry of branching in both wild-type and mutant kidneys is linked to the relative positioning of the tips, which is consistent with inter-tip suppression of bifurcation. Taken together, these findings show that elaboration of the ureteric tree occurs by way of a previously unrecognized but highly reproducible pattern of branching that contributes to the form and structure of the kidney.

We imaged ureteric trees from embryonic mice using optical projection tomography and quantified branching using Tree Surveyor software (Short et al., 2013). We analysed wild-type patterning in 32 C57BL/6J embryonic ureteric trees from E11.5 to E15.5. In observing the ureteric tree pattern, we noticed that each bifurcation led to significantly more tips from one child branch than the other, in a ratio that was broadly invariant. To quantify this, we defined the weight of a bifurcating branch point as the number of tips descending from it and defined its balance as the ratio between the larger and smaller child weights (Fig. 1A). We have previously proposed that the initial branches that establish the basic six-clade structure of the organ are distinct from further branching within clades (Fig. 1B) (Short et al., 2014). To examine this, we used our measures of weight and balance to study the root node as well as root branch points that establish the five ancestor and six-clade bifurcations that define the basic structure of the kidney. We refer to these as the 11 primary branch points and consider them separately from the other intra-clade branches that make up the vast majority of the ureteric tree (Fig. 1B). Branching at the ureter root node is significantly more symmetrical than in later generations (balance=1.18±0.18, mean±s.e.m.), with no anterior or posterior bias (P=0.495, paired t-test on log weight of anterior and posterior), whereas the next branches are substantially more asymmetrical than later generations (balance=2.21±1.20), principally because they define the clade structure of the organ (Fig. 1C). The remaining primary branch points have very similar balance (1.49±0.46, Fig. 1D). With the exception of the earliest developmental stage, balance at each of the primary branch points remains the same across developmental time (P>0.24 for linear regression against age as continuous or nominal variable, controlling for clade).

We next examined the balance of branches within individual clades, which make up the majority of the mature ureteric tree. In contrast to the variation we saw between many of the primary branches, the balance within clades was remarkably similar when comparing clades from all organs as a single dataset (balance=1.25±0.17, Fig. 1E): i.e. on average, the balance of any branch point within any clade is no different to any other branch point in another clade (P=0.98, ANOVA). The extensive reorganisation and branch resorption associated with the formation of the renal pelvis (Short et al., 2013) precludes the production of a ‘normally’ arborized ureteric after E15.5 (i.e. early branch generations are thereafter missing). However, to examine whether balance is maintained later in development, we investigated peripheral branching in four E16.5 kidneys. These organs have, on average, 1292 tips, which means that they are approximately one generation away from the cessation of branching. Using these data, we quantified balance in the three most distal branching generations (inverse generation 2,3,4 by our schema) and found no significant difference compared with younger organs (mean balance=1.252, P=0.1867 for a t-test on kidney means). This indicates that branching within clades occurs through the establishment of a constrained, asymmetrical and self-similar pattern of elaboration at least until the very latest stage of branching morphogenesis. The conservation of balance over time, irrespective of location within the organ, suggests that the elaboration of the ureteric tree occurs in a highly controlled manner after the initial establishment of basic organ pattern (through the primary branch points).

The conserved branching asymmetry we observe suggests that a consistent, locally operative mechanism defines a consistent pattern of branching morphogenesis that drives renal development. Although other biological systems such as the mammary gland have been successfully modelled using stochastic approaches (Scheele et al., 2017), after consideration these were found to not be appropriate in the case of the kidney, for reasons outlined later in the Discussion. Instead, we developed theoretical ‘tip state’ models to represent the evolution of the branching tree over time, with a view to describing our biological data. These models posit that a given tip exists in a state governed by rules that define whether it either matures or bifurcates as development progresses and what state it next assumes. The simplest such model is the unique non-trivial ‘single state’ model giving a family of trees with uniform bifurcation, which we call perfect. Because by definition this cannot model the systematic asymmetry observed in the analysis of our experimental models, we then considered ‘two-’ and ‘three-state’ models. There are two distinct two-state models and 34 distinct three-state models (see supplementary Materials and Methods for details and a classification of all models with up to three states).

We assessed all of these and found only three that have a single bifurcating state with non-terminating child states. Such generalized delay models have the desired property that balance converges to a single asymptotic value as weight increases. Two of these were considered plausible candidates for modelling the ureteric tree: a two-state model, which we term the ‘Fibonacci’, in which the weight of all trees and sub-trees are Fibonacci numbers (Fig. 2A); and a three-state model that we term the ‘half delay’, the behaviour of which is outlined in Fig. 2B. The Fibonacci model has asymptotic balance ≈1.618 (the golden ratio), while the half-delay model has asymptotic balance ≈1.325. The third delay model has asymptotic balance ≈2.148, which is outside the observed range in ureteric trees. In all cases, the smaller trees produced by these models have balance values that differ from the asymptote and from each other, owing to a discrete weight effect.

To define how these models of branching describe the behaviour of the ureteric bud in vivo, we compared the Fibonacci and half-delay models with our biological dataset using an optimal overlay algorithm (Lamberton et al., 2015), using the perfect model as a baseline. This approach computes a matching score between 0 and 1 for each pair of tree structures (Fig. 3A,B). We applied this method by clade, such that for every ureteric tree each of the six clades was overlaid with the closest tree by size from the selected model. We then calculated the overall matching score for the whole tree. When such a comparison was performed using our wild-type datasets, the best observed fit was the half-delay model that had a matching score between 0.74 and 1.0 depending on developmental stage; Fig. 3C shows an example of a subtree with an exact matching to the half-delay model. To understand the nature of the small amount of discordance between the model and in vivo branching of the ureteric bud, we undertook a program of pairwise matching between the model and wild-type data. For each developmental stage and clade type, we isolated the corresponding set of clades then identified the half-delay tree with best average matching score and added this tree to the set. We asked which member of this augmented set was best matched to all of the other members and then termed this the centroid (Lamberton et al., 2015). Using this approach we found that of the 36 cases (six clades and six developmental stages), the half-delay model was the centroid in 29 (80.6%), was the second ranked in six cases and third in only one (Fig. S1). This indicates that the model trees are ‘in the middle’ of the experimental ones that vary around them and that the half-delay model almost perfectly represents the majority of branching morphogenesis in the ureteric tree.

To define the normal developmental program that gives rise to the mature ureteric tree, we took advantage of our datasets to capture snapshots of the average branch pattern. The frequency of each subtree type was calculated across the entire dataset and the centroid for each weight identified as the single most representative structure (Fig. 4, top row). Although multiple tree structures occur at all weights above three, in each case the distribution is dominated by either one or two patterns (Fig. 4, for weights up to 12). In the latter instance, this always represents only a transient divergence. Strikingly, for the subtree weights that can be described by a half-delay tree model, the half-delay model dominates the distribution for that weight (see Fig. 4 weights 3, 4, 5, 7, 9 and 12). Similar calculations at weights 16 and 21 (the other half-delay examples) also found that that the centroid trees are the half-delay models of those weights (Fig. S2). This indicates a definable sequence of growth in the kidney, with the addition of a bifurcation at each step (green in Fig. 4) that is centred on the half-delay model.

One explanation for the asymmetry we have observed in whole clade and subtree analyses is that bifurcation of a given tip may be locally inhibited by its neighbours. If this were the case, we would expect the most closely related tips in the tree structure to provide the greatest local influence on growth and bifurcation, and hence expect a relationship between sibling node imbalance and the position of the next most closely related tips. To test this hypothesis, we analysed all pairs of sibling nodes with unequal weight (32 kidneys, 1963 pairs) and identified the common parent node p and its sibling, p′ (Fig. 5A,B). In the original pair of sibling nodes, we call the node that is spatially closest to p′ the inner sibling (b), which corresponds to the angle θ₁ being greater than θ₂ (Fig. 5A,B). If inter-tip suppression is driving asymmetry, we would expect the inner sibling to face stronger inhibitory signals and to have a lower weight (i.e. less branching). Analysis of these sibling pairs described such a scenario in 69.9% of cases (n=1963), supporting this hypothesis. In addition, nodes are not generally co-planar, as depicted in Fig. 5, meaning that a and b could be almost equidistant from p′ (Fig. 5B). To correct for this we defined near-planar branches where the inner sibling is more strongly ‘between’ p′ and a as those where the θ₁/θ₂ difference is above 45°. In these cases, the proportion in which the inner sibling weight is lower rises to 83.5% (n=571), a relationship that is apparent at all of the developmental stages we examined (>66% at each stage).

To test whether the delayed bifurcation of the inner sibling is just the result of generally slowed growth, we compared the branch lengths in pairs of near-planar sibling nodes (in Fig. 5A,B we compare the lengths from p to a and p to b). In cases where both siblings have bifurcated, the inner sibling is longer in 65.5% of cases (median difference 14.6 µm or 19.5%), indicating that a failure to branch resulted in continued lengthening. To further clarify this observation, we considered node pairs in which neither sibling has branched (i.e. a and b in Fig. 5A are both tips). Doing so allowed us to consider tip extension in isolation from bifurcation and in such cases little difference in length was seen: the inner child is instead shorter in close to half of cases (56.9%) and the median difference is much smaller than in siblings that have bifurcated (only 3 µm or 5%). Based on these findings, we conclude that the local cues that we propose to inhibit tip bifurcation have little or no impact on the rate of tip extension.

We next wished to address how the pattern of branching morphogenesis we observed during kidney development is established, with a view to understanding its genetic basis and robustness in the face of perturbation. On the basis that local cues are likely to influence patterning, we analysed mutants in which the environment of tips on the organ surface are altered. The first selected was bone morphogenetic protein 7 (Bmp7^−/−), an inhibitor of nephron progenitor cell differentiation (Tomita et al., 2013; Muthukrishnan et al., 2015) produced by cells in the UB epithelia and associated mesenchyme, and whose germline inactivation leads to renal hypoplasia (Dudley et al., 1995). Knockdown of expression in cultured organs results in fusion of ureteric tips and hence the molecule has been proposed to be essential in a tip-avoidance mechanism in the developing organ (Davies et al., 2014). The second was sprouty 1 (Spry1^−/−), a receptor tyrosine kinase inhibitor expressed in the ureteric bud tips (Zhang et al., 2001) that plays a central role in modulating cell signalling downstream of the Ret, a co-receptor for GDNF (Basson et al., 2005). Gene knockout results in precocious UB development and accelerated branching morphogenesis (Basson et al., 2005). Finally, we studied mice carrying heterozygous loss-of-function mutations in transforming growth factor β 2 (Tgfb2^+/−), a molecule expressed predominantly in the ureteric epithelium (Pelton et al., 1991). Tgfb2 haploinsufficiency results in hypoplasia characterized by reduced branching and developmental delay (Short et al., 2010).

To examine the relationship between tip spacing/environment and developmental patterning of the branching ureteric tree, we collected foetal kidneys from Bmp7^−/− (E14.5), Spry1^−/− (E13.5) and Tgfb2^+/− (E14.5) embryos and subjected them to OPT imaging and Tree Surveyor analysis as per the wild-type sample cohort used in our analysis of branch patterning (see supplementary Materials and Methods for basic metrics). Mutant kidneys of all three genotypes were found to be smaller than controls (Fig. 6A-C, Movie 1). To examine how this related to the physical environment in which the tips reside, we computed a surface area for each ureteric tree using a convex hull around the organ using the tip coordinates to establish the tree boundary (see supplementary Materials and Methods), which confirmed the hypoplasty in these organs (Fig. 6E). When tip environment was expressed as a surface area per tip ‘occupancy’, we observed no change in Tgfb2 kidneys (12,059±1453 µm² versus 10,396±351 µm², P=0.56), an increase in Bmp7 kidney (27,292±2844 µm² versus 11,783±489 µm², P=6.8×10⁻⁶) and a decrease in Spry1 kidney (6761±1038 µm² versus 12,192±445 µm², P=0.0044) (Fig. 6F). These observations are consistent with the respective hypoplasty and increased branching previously reported for these genotypes (Dudley et al., 1995; Basson et al., 2005). Importantly, differences in occupancy were still observed when mutant kidneys were instead compared with wild-type organs with comparable tip numbers (Fig. 6A-C, Movie 1), suggesting that such changes are not simple products of developmental delay but instead represent genetically encoded mediators of tree morphology (Fig. S3A). Notably, we observed no evidence for tip fusion or ‘collision’ in Bmp7 mutant kidneys, as has been reported in cultured organs (Davies et al., 2014). In fact, the surface area occupied per tip in the Bmp7 samples was considerably greater than in wild-type samples (see Fig. 6). Consistent with this observation, the average distance between tips was increased by 25.17 µm when compared with age-matched controls [95% CI (18.36,32.08), P=6.8×10⁻⁶] and by 17.42 µm when compared with wild-type kidneys with equivalent tip numbers [95%CI (9.94,24.84), P=2.0×10⁻⁴]. Differences in branching morphogenesis between cultured kidneys and organs dissected from mutant embryos are a feature we have noted before in several other mutants (Short et al., 2013). These may reflect artefactual products of anlage culture systems derived from tissue compression that interrupt the normal spatial arrangement of tissues and cells in the developing kidney.

Although the average distance is greater in Bmp7, we also considered how uniformly distributed the tips on the surface of each of the kidneys are for all genotypes. To determine this, we analysed the relative distance ratio of the closest and furthest neighbour tips from pairs of sibling tips. While this measure was unchanged when comparing Tgfb2 kidneys and controls (1.145684 versus 1.131798, P=0.1111, Fig. 6G), we noted significant differences between Spry1 and Bmp7, and age-matched controls (1.161907 versus 1.120212, P=0.0485 and 1.218161 versus 1.122488, P=1.361×10⁻⁵, Fig. 6H,I). We interpret this to mean that the distribution of tips on the surface of these organs is less even.

Having identified mutants with altered tip spacing and distribution, we employed overlay metrics to determine whether these changes in the tip environment and in branching dynamics resulted from, or were the consequence of, perturbation of the half-delay patterning we had described in wild-type organs. The average balance was calculated for each kidney (Fig. S4) and each mutant group was compared with the age-matched wild-type control group (Fig. 7A). The Tgfb2 mutants show no detectable change in balance, consistent with our previous overlay analysis (Lamberton et al., 2015), but increases in branching asymmetry were observed in the Bmp7 and Spry1 mutants. Comparisons were repeated using size-matched controls (by tip number), with similar results (Fig. S3B). We then repeated the centroid analysis previously employed to characterise the wild-type dataset, examining the Bmp7 and Spry1 mutants for each of the six clades. For each mutant group, the selected clade was isolated from each kidney, and the resultant set of clades was augmented with the closest matching model tree (in this case either half delay or Fibonacci) before finding the centroid. In both cases, the half-delay model was the centroid in only 33% of clades (two out of six) as opposed to >80% of wild-type organs at the same developmental stage (29 of 36 cases, Fig. S1). In contrast, when the Fibonacci model was used to classify the clades, it was the centroid in 50% of cases for Spry1 (three out of six) and ∼83% of cases for Bmp7 (five out of six). This modification in renal patterning is also reflected in the distribution of branch balances for different sub-tree sizes (Fig. 7B), where a greater imbalance in the mutant trees is evident. Hence, in both the Spry1 and Bmp7 kidneys, where there is a significant change in the spacing of tips on the organ surface, there is a definite shift from the wild-type half-delay model towards the Fibonacci model. This itself represents a change in the relationship between the extension and bifurcation rates in the mutant trees.

Branching morphogenesis is a commonly employed developmental mechanism that establishes the airways of the lungs, the ductal network of glandular organs and the collecting duct network of the kidney. In all cases, the epithelial networks present in these organs are established by way of the arborisation of an initial epithelial bud, usually within a supporting population of mesenchymal cells. The factors that mediate this elaboration and the manner in which it occurs are only partially defined. In the lung, it has been proposed that a rigid and genetically hardwired program of three different branch types combine to form the airways. However, this pattern of growth appears to differ fundamentally from the kidney, as there is no recognizable and reproducible structure of the branched ureteric tree (Short et al., 2014). Our previous studies have identified similarity between different kidneys in the form of a higher order clade structure and also in the extent of branching (Short et al., 2014). In this study, we have sought to understand whether (and to what extent) this reflects a deeper underlying patterning that governs the development of the organ.

We have established a mathematical growth model of ureteric tree branching in which a stereotypical early structure is succeeded by a branching pattern that is locally similar regardless of the position of a branch within the ureteric tree and across developmental stages. This pattern is characterized by a conserved level of branching asymmetry, suggesting a common growth process acting at a local level. This is consistent with our previous quantitative analysis of branch sizes and geometry (Short et al., 2014), which similarly shows substantial underlying regularity, with minor local variation. The family of generative tip-state models that we have developed comprise tip populations that either bifurcate or mature solely according to their current state. These can be conceived more generally as delay models (see supplementary Materials and Methods), in which there is a characteristic difference in maturity or time to bifurcation between the two child nodes formed by each bifurcation. Consequently, they possess the feature that the balance at each branch point will converge to a fixed asymptotic value as growth continues, and are well suited to modelling branching under a common locally operative growth process. We have shown that a three-state delay model (which we term ‘half delay’) recapitulates, with a very high degree of accuracy, the pattern of branching established by imaging biological samples. This encompasses measures of balance that are shared between and within the clades of the organs. Importantly this model accurately describes the structure, both on the global (Fig. 3) and local (Fig. 4) scale.

Our analysis suggests that an initial idiosyncratic phase of stereotypic UB branching occurs early in development and establishes the clades of the organs. The organization and balance of these branches is clearly distinct from later bifurcation events. At this point, the UB invades a discrete and relatively large and uniform region of SIX2⁺ GDNF-producing mesenchyme. We propose that the early branching pattern and the establishment of clades occurs in response to the pre-existing shape and relative uniformity of this field of cells. However, by E12 the SIX2⁺ population begins to condense around the growing tips to form discrete cap mesenchyme populations, which are then a feature of organ development until the termination of nephrogenesis at P4 in the mouse. During this period of development, the tips are confined to, and evenly arranged within, a thin zone on the surface of the growing kidney. We propose that this environment is conducive to the establishment of locally operant cell signalling systems that dictate the asymmetric pattern of branching morphogenesis we observe in our analysis. Although the surface area and tip number increase exponentially, the spacing between tips declines only modestly in synchrony with the reduction in the size of the tip and cap niches (Short et al., 2014). The fact that tip spacing is relatively constant indicates coordination between tip extension and branching: the two primary ureteric tree growth processes. A natural hypothesis is that bifurcation of the growing tips is suppressed to some extent by neighbouring tips and we propose that such a mechanism further contributes to the patterned branching we observe. The hypothesized mutual suppression of bifurcation by the tips is supported by our geometrical analysis. Bifurcation is delayed (but growth continues) in tips located between the most closely related branches in the tree structure (node b in Fig. 5), which we expect to have, on average, more and closer neighbouring tips than their siblings (node a in Fig. 5). In this way relatedness in the tree structure provides a useful proxy for tip spacing. Furthermore, this comparison of sibling branches demonstrates that the mutual suppression of bifurcation by tips may explain the systematic and conserved branching asymmetry that marks the branching structure.

Recent studies have demonstrated that the hierarchical patterning of the mouse mammary gland can be modelled with a locally operative, time invariant, stochastic growth process, in which growing tips bifurcate or cease growth with near-equal probability (Scheele et al., 2017). This report further suggested that a similar process may drive renal branching morphogenesis. We find little support for such a proposal. For example, there is no evidence of stochastic termination of bifurcation or extension – branches do not terminate ‘within’ the kidney and tips are all, by and large, bifurcating on the organ surface at any given time (see Fig. 6 and Movie 1). Furthermore, cell proliferation in all of the tips of the branching UB is actively maintained until the cessation of nephrogenesis at postnatal day 4 (Short et al., 2014), in a manner that relies on signals from the cap mesenchyme (Cebrian et al., 2014). The possibility that the observed asymmetry may arise from a random process was also considered. The clearest evidence against this hypothesis (and on derived models in which the timing of random branching is shaped by a specified probability distribution) is found in the geometrical analysis of sibling subtree clusters (see Figs 3 and 4). As we can predict which sibling subtree is most likely to be larger based on their positions relative to related tips, a purely random process can be excluded. However, ultimately the strongest argument for this approach is that we were able to provide an excellent fit to the centre of the distribution of real data (see centroid analysis). In these analyses the half-delay model trees were generally a better fit to the experimental data than were any of the individual experimentally determined trees. The remaining option for a stochastic model would be to add some degree of randomness to the half-delay model. Although this might be an accurate representation of the observed biological variation around the identified stereotypical pattern, the lack of any systemic difference between our deterministic model and the data indicated that fitting such a model would not further improve the description.

Taken together, these observations are entirely inconsistent with a model based on stochastic cessation of branch growth and/or bifurcation. Instead, we propose that there are much broader similarities between branching morphogenesis in the kidney and lung. Unlike the mammary gland, branching morphogenesis in these organs is particularly characterized by the involvement of highly specialized mesenchymal cell populations that produce factors that drive cell proliferation and branching of their respective epithelial networks. In support of this notion, we have found that the tip state models developed in this study are able to generate hierarchical structures like the domain branching seen in lungs (see supplementary Materials and Methods: three state models). Although further modelling and comparison with biological samples will be required to determine the full extent to which such models can describe pulmonary development, our findings suggest a commonality between branching in these organs and highlight considerable differences with the mammary gland.

Our analyses of mutant kidneys say much about the nature of patterning that governs the normal development of the kidney. They imply that alterations in the rate of organ growth do not necessarily perturb patterning. While haplo-insufficiency for Tgfb2 significantly abrogates kidney growth, by ∼24 h (Lamberton et al., 2015), we observe no differences in the pattern of systematic asymmetry in branching identified in this study. In most respects, these organs, although small, are perfectly formed – suggesting that genetic control of organ growth is not intrinsically linked to patterning of the branching UB. In contrast, changes in Bmp7 and Spry1 mutants result in changes in branch patterning from a half delay to a pattern more closely approximating a Fibonacci model. Although the surface areas per tip of each model are larger (Bmp7) and smaller (Spry1), both models share the feature of more uneven tip distribution on the organ surface. This suggests that the normal factors governing tip distribution are somehow ameliorated in both organs. This might be through a reduction in suppressive actions of a protein that prevents the incursion of one tip into the field of another (as may be the case for Bmp7) or in the loss of factors that normally control the unchecked bifurcation of tips (as may be the case for Spry1). In both cases, a change in the uniformity of tip spacing is associated with altered branch patterning. At this point, we cannot say with certainty whether the alteration in patterning is a cause or a consequence of the change in the uniformity of tip distribution. However, given that the patterning is effectively determined by the decision of a tip to bifurcate or grow (which is likely locally determined), we favour the former model in which the uniform distribution of tips on the organ surface dictates the systematic asymmetry we observe in the ureteric tree.

What biological factors might contribute to the half-delay models that describe patterning of the elaborating ureteric tree? The different ‘nodes’ in this scheme effectively represent points at which a given branch bifurcates or grows. Such states effectively represent different degrees of developmental delay (see supplementary Materials and Methods for approaches to modelling different delay to bifurcation time ratios). It is likely that both biophysical and molecular cues contribute to the decision of a given tip to bifurcate or to extend. The size of the organ surface effectively establishes the field in which such tip-tip interactions can occur and although it might be expected that this is defined by the extent of ureteric branching, our findings demonstrate that this is not the case: increased and decreased branching do not necessarily imply increased or decreased organ size. At the level of the tip, the bounding basement membrane likely constitutes a structure that is able to sterically impact adjacent tip environments. It is tempting to attribute similar features to the surrounding cap mesenchyme that (at least in section) appears as a ramified buffer between neighbouring tips. However, recent studies suggest that this cell population is dynamic and motile, which may reduce its steric impact on neighbouring tips (Combes et al., 2016). Equally, the role of the stromal cells – which sit between cap/tip populations – has been relatively poorly studied, although they are known to regulate branching morphogenesis in a non-cell-autonomous manner (Paroly et al., 2013; Mukherjee et al., 2017). Studies to assess the in vivo forces and tensions created by juxtaposed cap/tip/stromal niches and their respective cells are required to assess their biophysical impact on patterning. Considerable recent conceptual and experimental advances have extended our understanding of the cellular events that underlie bifurcation in branching epithelia (Ihermann-Hella et al., 2014; Menshykau et al., 2014; Huebner et al., 2016; Lin et al., 2017). Linking these to the emerging picture of shared and organ-specific growth factor production in the tips of branching epithelial networks (Rutledge et al., 2017) will therefore be important for integrating how these combine to direct the exquisite patterning we describe in the developing kidney at a whole organ level.

Embryonic mice from heterozygous Tgfb2^tm1Doe, homozygous Bmp7^tm1.2Dgra, homozygous Spry1^tm1Jdli and wild-type C57BL/6J kidneys were collected, stained and imaged as previously described (Combes et al., 2014) with ages and sample numbers outlined in Table S1. Limb staging was used to broadly confirm age in most of the embryos and also to subclassify them based upon developmental progression (Wanek et al., 1989). Kidneys were stained using anti pan-cytokeratin (Abcam, AB11213; 1:300) and anti Trop2 (R&D Systems, AF1122; 1:100) with Alexa Fluor secondary antibodies (Life Technologies), followed by OPT imaging (Bioptonics). All animal procedures complied with standards set under Australian, US and Canadian federal and state guidelines for animal welfare, and experiments were approved by the relevant Animal Ethics Panels of Monash University, Columbia University and the University of Alberta.

Reconstructed ureteric tree datasets were mapped and quantified in 3D using Tree Surveyor software (Short et al., 2013). Tree Surveyor software is available to academic and not-for-profit researchers. GraphML data was exported, providing a representation of each tree in the form of nodes (branch points and tips) and edges between them (branch segments), annotated with information on size and position. The root node of the graph is the point where the ureter first branches, and defines an orientation of each edge from a start node, which is closer to the root, to an end node, which is further away. The start node is the parent of the end node, and the end node is the child of the start node. Each of the 11 primary branch points (Fig. 1B) were annotated on the ureteric tree graphs based on their relative position in the kidney the graph was derived from, also identifying the six-clade structure previously described (Short et al., 2014).

Graph overlays and pairwise comparisons were performed as previously described (Lamberton et al., 2015) using MATLAB and R code. All other analyses described were performed in R, including generation of tip-state model tree structures used in overlay analysis with real data. See supplementary Materials and Methods for details of tip state model definitions, properties and enumeration. Example model trees used in overlay analysis were generated in R.

We acknowledge the contributions of Frank Costantini and Jonathan Licht in providing the Sprouty 1 mouse embryos. We thank the Monash Micro Imaging and Monash Animal Research Platforms for assistance with imaging and animal husbandry, respectively. The authors also thank Professor Melissa Little and Dr Alexander Combes for interesting and useful discussions during the preparation of this paper.