1 
 
 
Identification of mammary stem cells and their dynamics 1 
during branching morphogenesis 2 
 3 
 4 
Authors: Colinda L. G. J. Scheele* 1, Edouard Hannezo* 2,3,4, Mauro J. Muraro1, Anoek 5 
Zomer1, Nathalia S. M. Langedijk1, Alexander van Oudenaarden1, Benjamin D. 6 
Simons2,3,4, Jacco van Rheenen1  7 
 8 
Affiliations: 9 
1Cancer Genomics Netherlands, Hubrecht Institute-KNAW & University Medical Centre 10 
Utrecht, Utrecht, the Netherlands 11 
2Cavendish Laboratory, Department of Physics, University of Cambridge, Cambridge, 12 
United Kingdom 13 
3The Wellcome Trust/Cancer Research UK Gurdon Institute, University of Cambridge, 14 
Cambridge, United Kingdom 15 
4The Wellcome Trust/Medical Research Council Stem Cell Institute, University of 16 
Cambridge, United Kingdom 17 
*  These authors contributed equally to this work 18 
 19 
Corresponding authors:  20 
JvR: j.vanrheenen@hubrecht.eu  21 
BDS: bds10@cam.ac.uk.  22 
2 
 
 
Abstract 23 
At birth, the mouse mammary gland consists of a rudimentary tree that, by ductal 24 
elongation and terminal end-bud (TEB) bifurcation, develops during puberty into a 25 
highly branched epithelial network. Based on morphology and expression of stem 26 
cell markers, mammary stem cells (MaSCs) are thought to reside in or near the 27 
TEBs. However, in the absence of exclusive stem cell markers, the true identity of 28 
MaSCs, including their location, multiplicity, dynamics and fate during branching 29 
morphogenesis, has yet to be defined. Here, we developed a multidisciplinary 30 
approach based on whole gland reconstruction, genetic lineage tracing, 31 
mathematical modelling and single cell sequencing to define the identity and fate 32 
of MaSCs during mammary epithelial morphogenesis. At the organ scale, 33 
morphogenesis happens in a time-invariant manner driven by TEB bifurcation and 34 
termination, occurring stochastically with near equal probability, leading to a 35 
heterogeneous network topology. We show that the majority of TEB cells function 36 
as highly proliferative, lineage committed, MaSCs that are heterogeneous in their 37 
expression profile and their short-term contribution to ductal extension. Yet, 38 
through cell rearrangements following TEB bifurcation and elongation, each MaSC 39 
within this heterogeneous population is able to contribute actively to the growth of 40 
the epithelial network. Together, our study shows that the behaviour of MaSCs is 41 
not directly linked to a single expression profile. Instead, morphogenesis of the 42 
mammary epithelium relies upon lineage-restricted heterogeneous MaSC 43 
populations that function long-term as single equipotent pools.  44 
 45 
3 
 
 
Main Text 46 
During puberty, ductal elongation is mediated by cell proliferation at TEBs (Fig. 1a and 47 
Extended Data Fig. 1)1-3. This proliferation would quickly subside without the support of a 48 
distal self-renewing cell population, suggesting that MaSCs must reside either close to or 49 
within TEBs4,5. MaSCs have been associated with various morphologically distinct cells6-50 
8, including cap cells that line TEBs1,2, but their true identity has yet to be defined. General 51 
basal markers, such as Keratin 14 (K14) and K8, have been used to study the luminal 52 
and basal lineages of MaSCs9-11, but could not be used to exclusively label MaSCs. 53 
Attempts to define specific markers, such as s-Ship, Axin2, Lgr5, and Procr, have led to 54 
ambiguous results whereby the genes were either not exclusively expressed in MaSCs 55 
or only expressed by a subset of MaSCs9-18. Importantly, as MaSCs are expected to 56 
contribute throughout development, their potency, hierarchical organization and fate must 57 
be assessed from functional assays that record their behaviour over the long term. 58 
Moreover, although much work has been devoted to the regulation of branching19, 20, it 59 
remains unclear how MaSCs specify the macroscopic architecture of the gland. Here, to 60 
address these questions, we have combined unbiased genetic lineage tracing and whole-61 
gland mapping, to identify the location and number of MaSCs, and to reconstruct 62 
quantitatively how the fate of MaSCs translates into organ structure.  63 
 64 
Non-stereotypical mammary branching pattern 65 
To define the development of the ductal epithelium at the organ scale, 4th and 5th 66 
mammary glands were isolated from mice that had reached adulthood (8 weeks) and 67 
stained for K14. To reconstruct the ductal tree, and resolve basal and luminal cell identity, 68 
4 
 
 
we obtained images at the organ-scale by stitching together >200 high-resolution Z stacks 69 
(50 steps with 5µm step size) (Fig. 1b). The resulting structure of the epithelial network 70 
was specified through “level indices”, defined as the number of branches between a given 71 
duct and the rudimentary tree (Extended Data Fig. 1d). In contrast to regularly branched 72 
epithelia such as the lung19, the mammary gland is characterised by morphological 73 
heterogeneity20, with both large and small subtrees (Fig. 1c). This raised the question of 74 
whether the inter-subtree heterogeneity derives from the early specification of 75 
heterogeneous TEBs with defined proliferative potential, or whether mechanical or other 76 
local cues21-23 condition disparate fates of an otherwise homogenous population of TEBs.  77 
To address this question, we looked for evidence in the growth characteristics of the 78 
reconstructed ductal trees. First, quantitative comparison of the branch length 79 
distributions between the ductal trees of 5 versus 8 week-old mice suggests that 80 
branching occurs nearly exclusively from TEB bifurcation, as opposed to side-branching 81 
from existing ducts (Extended Data Fig. 1e and Supplementary Theory). EdU labelling at 82 
various stages of puberty confirmed proliferation at the TEBs and not within existing 83 
ducts, as would be expected for side-branching (Extended Data Fig. 1c). Interestingly, we 84 
also noted that proliferation did not subside gradually throughout the gland (Extended 85 
Data Fig. 1f). Instead, during the course of development, an increasing number of TEBs 86 
become completely non-proliferative, while the remainder maintain their full proliferative 87 
activity until, by the end of puberty, all TEBs have exited cycle (Extended Data Fig. 1g). 88 
This implies that MaSCs within each TEB act collectively to regulate either TEB 89 
bifurcation or termination. Indeed, analysis of the frequency of these events from the 90 
reconstructions shows that both occur randomly with nearly equal probability at any point 91 
5 
 
 
along the ductal tree (Fig. 1d and Extended Data Fig. 1h). Importantly, using only the 92 
average balance between TEB bifurcation and termination, we could predict quantitatively 93 
the subtree persistence (Fig. 1e, left panel), the subtree size heterogeneity (Fig. 1e, right 94 
panel), as well as the diminishing fraction of proliferative TEBs over time (Extended Data 95 
Fig. 1g and Supplementary Theory). Together, these findings point at a simple design 96 
principle in which all TEBs have the same growth potential and follow a time-invariant 97 
stochastic rule throughout pubertal mammary development. Consistent with this 98 
behaviour, we found that the average length and width of ducts, as well as the number of 99 
proliferative cells per TEB, remain level-independent throughout development (Extended 100 
Data Fig. 1i,j).  101 
 102 
Lineage tracing reveals the location and number of MaSCs 103 
Although the organ-level analysis provides insight into the collective behaviour of 104 
MaSCs within a given TEB, it does not reveal the potency, multiplicity and fate behaviour 105 
of individual MaSCs. To characterize the behaviour of individual MaSCs, we developed 106 
an unbiased clonal lineage tracing approach based on a ubiquitous promoter using R26-107 
CreERT2;R26-Confetti mice, allowing the fate of cells to be traced at single cell 108 
resolution24 (Extended Data Fig. 2a). Upon injection of an ultra-low dose of Tamoxifen 109 
(0.2mg/25g body weight, Extended Data Fig. 2b) at the onset of puberty, at a level that 110 
does not interfere with morphogenesis (Extended Data Fig. 3 and Supplementary 111 
Theory)10, 25, Cre-mediated recombination of the confetti colour randomizer was 112 
stochastically induced in a minority of cells (<1%) present in the rudimentary ducts, and 113 
mice were sacrificed either mid-puberty (5 weeks) or on reaching adulthood (8 weeks). 114 
6 
 
 
Whole mount reconstruction of the mammary gland revealed isolated confetti-labelled 115 
cells from all four confetti colours in the ducts of the rudimentary tree, consistent with the 116 
induction of non-proliferative cells (Extended Data Fig. 4a,b). At higher levels of the 117 
epithelial network, which were formed after induction, we found that subtrees were either 118 
entirely unlabelled, or consisted of a small fraction of labelled cells bearing only a single 119 
confetti colour from either the luminal or basal lineage, but not both (Fig. 2a,b and 120 
Extended Data Fig. 4c,d). Since the labelling density was extremely low (only 1 in 3 glands 121 
were labelled, n = 16 glands), we reasoned that each clone must have originated from a 122 
single founder cell. Moreover, since these cells contribute during the entire course of 123 
pubertal ductal development (Fig. 2c), we concluded that the founder cell must belong to 124 
the self-renewing MaSC pool. Importantly, longer-term tracings showed that clones 125 
derived from the same assay could contribute to the various stages of remodelling in 126 
adulthood, including the oestrous cycle-driven formation of small lobuloalveolar structures 127 
and pregnancy-driven formation of secretory alveoli (Extended Data Fig. 5a-d).  128 
To quantify the dynamics of clonally labelled cells, we induced a red-to-green switch 129 
at clonal density in the mammary glands of R26-CreERT2;R26-mTmG mice at the onset 130 
of puberty. Next, after two weeks of tracing, we imaged these mammary glands at high 131 
resolution for 150 minutes using multiphoton intravital microscopy through a mammary 132 
imaging window. In line with EdU labelling (Fig. 1a and Extended Data Fig. 1), intravital 133 
imaging revealed that, in contrast to the dynamic and proliferative TEB cells, ductal cells 134 
do not actively proliferate or migrate (n = 3 mice, 17 TEBs, 12 ducts) (Fig. 2d and 135 
Supplementary Video). Consistent with this observation, the size of labelled ductal cell 136 
clusters deposited between 3 and 5 weeks was identical to that found between 3 and 8 137 
7 
 
 
weeks (Extended Data Fig. 2c-f). Together, these findings imply that ductal cells provide 138 
a “frozen” record of MaSC fate once they are deposited, and that spatial information can 139 
be translated into temporal information on the past behaviour of MaSCs.  140 
Using this strategy, we sought define the location of MaSCs. In mid-puberty glands, 141 
we scored the position of single confetti+ cells localized either in the TEBs or the 142 
immediate ductal region (Fig. 2e, top panel and Extended Data Fig. 5e-h). We found that 143 
all singly isolated confetti+ cells were localized within the TEB (Fig. 2e, bottom panel) and, 144 
since these cells had given rise to confetti+ cells in previously formed ducts, we concluded 145 
that these must be MaSCs. To determine the fraction of cells in the TEB that were MaSCs, 146 
we analysed subtrees in clonally labelled glands. At clonal density labelling, we reasoned 147 
that, in the TEBs where tracing was initiated, only one out of the total number, N, of 148 
MaSCs was confetti labelled so that, in the ductal tree that was grown afterwards, only a 149 
fraction of ductal cells (1/N on average) should become labelled with a confetti colour 150 
(Extended Data Fig. 6a). By analysing this fraction based on K14 and E-cadherin staining 151 
of whole mount mammary glands (Extended Data Fig. 6b,c), we inferred an estimate of 152 
170±30 (mean±s.e.m.) luminal MaSCs and 90±30 (mean±s.e.m.) basal MaSCs per TEB 153 
(Fig. 2f). Since these values equate, respectively, to the approximate total number of 154 
luminal and basal cells present in TEBs (Extended Data Fig. 6d,e), this suggests that 155 
most TEB cells function as basally or luminally committed MaSCs contributing to the 156 
entire pubertal development of the ductal tree.  157 
 158 
Molecular characterization of MaSCs using single cell mRNA sequencing 159 
8 
 
 
Having identified functionally the location and number of MaSCs, we next aimed 160 
to characterize their molecular signature by single cell mRNA sequencing. Individual 161 
ducts and proliferative TEBs were micro-dissected from 5 week-old MMTV-Cre;R26-loxP-162 
stop-loxP-YFP mice and single mammary epithelial cells were sorted by isolating Lin-163 
(CD45-, CD31-, CD140a-) YFP+ cells using flow cytometry. Since epithelial cells may not 164 
be labelled equally in this model, a slight bias toward high YFP expression could not be 165 
avoided. Using the StemID algorithm26, we identified computationally 9 clusters of cells 166 
(n = 91 cells in total, Fig. 2g and Extended Data Fig. 7a). The absence of correlation 167 
between cell cycle genes and clusters showed that they were not associated with distinct 168 
cell cycle states (Extended Data Fig. 7b). Based on expression of canonical basal cells 169 
markers such as K5, K14 and smooth muscle actin, and luminal markers, such as K8, 170 
K18, and K19, we could assign clusters to the luminal or basal lineage (Fig. 2g and 171 
Extended Data Fig. 7c,d). Previously described MaSC markers, such as Axin2 and Lgr5, 172 
were lowly expressed and detected randomly in both ductal and TEB cells in multiple 173 
clusters (Extended Data Fig. 7e). Surprisingly, despite our finding that the majority of 174 
basal and luminal TEB cells function as MaSCs, differential gene expression profiling of 175 
the clusters suggested that both luminal and basal MaSCs can vary in their transcriptional 176 
program (Fig. 2g and Extended Data Fig. 7f,g). For example, Csn3 and Cst3, which are 177 
related to milk production27, are more highly expressed in luminal clusters 2 and 3 when 178 
compared to luminal cluster 1 (Extended Data Fig. 7f). Similarly, Myh11 and Mylk, which 179 
are involved in contractile activity28, are more highly expressed in basal clusters 6 and 9 180 
when compared to basal cluster 7 and 8 (Extended Data Fig. 7g). Together, these 181 
9 
 
 
observations question whether the various transcriptional programs reflect a spectrum of 182 
primed MaSC states within the TEB. 183 
 184 
MaSCs are functionally heterogeneous in the short-term 185 
To test whether molecular heterogeneity of MaSCs might translate to functional 186 
heterogeneity, we analysed in mid-puberty glands the short-term capacity of single 187 
confetti+ MaSCs to contribute to ductal growth by analysing the size of the subclone in 188 
the adjacent duct (Fig. 3a). Interestingly, we found that only single confetti+ MaSCs 189 
localized at the border of the TEB (but not at the tip) contributed to the growth of the 190 
adjacent duct (Fig. 3b). Consistently, the number of confetti+ cells in the adjacent duct 191 
was correlated with the number of similarly labelled confetti+ MaSCs at the border of the 192 
TEB, but not at the tip (Fig. 3c). Notably, this spatial advantage of border over tip MaSCs 193 
to contribute to ductal expansion is not caused by differences in proliferative activity since 194 
the EdU labelling intensity following a 3 day chase is diluted by cell division to a similar 195 
degree in both regions (Fig. 3d,e). These results suggest an organization of MaSCs in 196 
the TEB that mirrors that of adult stem cells in the intestinal crypt29. In both cases, cells 197 
positioned at the border of the “niche” experience a short-term bias towards 198 
differentiation, while cells at the base or tip are primed for renewal30. From this short-term 199 
positional bias, we expect “intermittency” in the ductal clonal record. Consistently, 200 
reconstruction of clones across multiple branch levels showed that MaSCs could cease 201 
contributing to ductal expansion over extended periods of time (Extended Data Fig. 4e 202 
and Supplementary Theory). However, we cannot rule out the alternative possibility that 203 
10 
 
 
these interruptions could be caused by a small population of cells (<5-10%) that remains 204 
quiescent in the TEB. 205 
 206 
MaSC mixing resolves short-term heterogeneity into homogenous long-term 207 
behaviour 208 
To determine whether the short-term bias of border cells to contribute to ductal 209 
expansion reflects long-term differences in MaSC potential, we analysed the evolution of 210 
subclone sizes (defined as the density of lineage labelled cells in a given branch) as a 211 
function of level index (which serves as proxy for time) (Fig. 4a). If fate bias would persist 212 
long-term, it is expected that just a few clones would grow large, and that the subclone 213 
size distribution would change over time. To our surprise, the distribution of subclone 214 
sizes, when rescaled by the average, was equivalent at all branch levels and consistently 215 
well-fit by an exponential (Fig. 4b and Extended Data Fig. 2g). Such “scaling” behaviour 216 
is a hallmark of neutral dynamics and shows that, long-term, MaSCs function as an 217 
equipotent pool (Supplementary Theory). 218 
To understand how functional short-term heterogeneity of MaSCs can be consistent 219 
with homogenous long-term behaviour, we used a modelling approach to predict the 220 
dynamics of MaSC clones during TEB bifurcation. Our analysis showed that, during the 221 
repetitive process of TEB bifurcation and ductal elongation, MaSCs are randomly 222 
segregated, leading to a drift of subclone size within TEBs until subclones become either 223 
lost or “fixed” (Fig. 4c and Extended Data Fig. 8a) at a rate dependent only on the degree 224 
of cell mixing between two consecutive bifurcations. To challenge the model, we tested 225 
whether it could predict the fate of multiple differentially coloured MaSCs following serial 226 
11 
 
 
rounds of branching. In 3 week-old mice, we induced the expression of confetti colours at 227 
a slightly higher density (1.5mg/25g Tamoxifen), that does not influence the branching 228 
morphogenesis (Extended Data Fig. 3), resulting in mosaic differential colour-labelling of 229 
two or three independent MaSCs per TEB. Consistent with the model, we found that ducts 230 
formed soon after labelling showed a mixture of different confetti colours, while ducts 231 
produced late in development became either mono-coloured or devoid of any confetti+ 232 
cells (Fig. 4d and Extended Data Fig. 8a and 9). Moreover, we found that the average 233 
number of labelled cells per duct in both basal and luminal clones increased linearly as a 234 
function of branch level, as predicted by a “well-mixed” model in which MaSCs are 235 
randomly segregated during bifurcation, independent of their relative position within the 236 
TEB (Fig. 3a, Extended Data Fig. 8b,c and Supplementary Theory). Crucially, we found 237 
that the labelled cell fraction remained approximately constant over time (i.e. branch level) 238 
post-induction (Extended Data Fig. 2c), confirming quantitatively the representativeness 239 
of labelling throughout pubertal growth. These features also held true in the 3 to 5 and 5 240 
to 8 weeks tracings (Extended Data Fig. 2e,f). Finally, we performed multiphoton intravital 241 
microscopy of TEB extension and confirmed that cells indeed mix and exchange position 242 
in TEBs during the growth phase (Fig. 2d and Supplementary Video). Together, these 243 
results explain how heterogeneous populations of lineage restricted basal and luminal 244 
MaSCs can behave long-term as homogeneous equipotent pools.  245 
Finally, in line with previous observations15, we noticed that the offspring of the MaSCs 246 
do not form a contiguous stream, but show a clustered pattern, i.e. uni-coloured confetti+ 247 
clusters appear alternately with clusters of non-labelled cells (Extended Data Fig. 8e). 248 
The size of the clusters from basal and luminal lineages were the same (Extended Data 249 
12 
 
 
Fig. 8f). This, in combination with our observation that the ratio of basal and luminal 250 
MaSCs reflects the ratio of basal and luminal cells deposited in the duct (Fig. 2e and 251 
Extended Data Fig. 3h), suggest that individual luminal and basal MaSCs have the same 252 
potential to produce offspring with limited proliferative capacity, resulting in cluster 253 
formation (Extended Data Fig. 8f).  254 
 255 
Discussion 256 
Based on these findings, we propose a model of mammary morphogenesis in which 257 
lineage-restricted pools of MaSCs remain at the traveling TEBs, leaving behind offspring 258 
(Fig. 5). Although only MaSCs positioned at the border of the TEB contribute transiently 259 
to ductal expansion, cell rearrangements during TEB bifurcation allow positional biases 260 
to be reassigned, resulting in long-term equipotency of the MaSC pools. Whether the 261 
observed transcriptional heterogeneity of TEB cells is a manifestation of evolving 262 
positional biases or a signature of discrete MaSC subpopulations in the luminal and basal 263 
lineages remains an interesting open question. Importantly, this study emphasizes that 264 
the behaviour and even the identity of stem cells cannot be directly linked to a single 265 
expression profile or specific markers, but instead must be defined functionally. 266 
Lastly, at the organ scale, our studies show that MaSCs act cooperatively to regulate 267 
their proliferative activity within individual TEBs, with TEB branching and differentiation 268 
occurring randomly with nearly equal probability leading to a heterogeneous network 269 
structure. Although the mechanisms that regulate such balance remain undefined, it is 270 
interesting to question whether other ductal tissues might conform to the same growth 271 
characteristic. Intriguingly, analysis of the E16 murine kidney showed that its complex 3D 272 
13 
 
 
ductal network topology could also be explained by the same model paradigm (Extended 273 
Data Fig. 10) suggesting that the underlying mechanism of branching morphogenesis 274 
may be conserved across different ductal tissues.  275 
14 
 
 
Methods summary 
Mice. All mice were from a mixed background, were housed under standard laboratory 
conditions, and received food and water ad libitum. All experiments were carried out in 
accordance with the Animal Welfare Committee of the Royal Netherlands Academy of 
Arts and Sciences, The Netherlands. R26-Confetti24;R26-CreERT2, R26-TdTomato;R26-
CreERT2, and R26-mTmG;R26-CreERT2 mice were injected intraperitoneally (IP) with 
Tamoxifen (Sigma Aldrich) diluted in sunflower oil to activate Cre recombinase at 3 weeks 
of age. To label only one cell per duct with a confetti colour, R26-Confetti mice were 
injected with 0.2mg/25g Tamoxifen. To label multiple stem cells per TEB, mice were 
injected with 1.5mg/25g at 3 weeks. The clonal dose for R26-mTmG and R26-TdTomato 
reporter mice was 0.2mg/25g and 0.05mg/25g Tamoxifen respectively. Lineage-traced 
mice were sacrificed at mid-puberty (5 weeks-old) or at the end of puberty (8 weeks-old) 
and all the mammary glands were collected. 
 
Whole mount imaging of mammary glands. Imaging of whole mount mammary glands 
was performed using a Leica TCS SP5 confocal microscope, equipped with a 405 nm 
laser, an argon laser, a DPSS 561 nm laser and a HeNe 633 nm laser. Different 
fluorophores were excited as follows: DAPI at 405 nm, CFP at 458 nm, GFP at 488 nm, 
YFP at 514 nm, RFP at 561 nm, and Alexa-647 at 633 nm. DAPI was collected at 440-
470 nm, CFP was collected at 470-485 nm, GFP was collected at 495-510 nm, YFP was 
collected at 540-570 nm, RFP was collected at 610-640 nm and Alexa-647 was collected 
at 650-700 nm. All images were acquired with a 20x (HCX IRAPO N.A. 0.70 WD 0.5 mm) 
15 
 
 
dry objective using a Z-step size of 5 µm (total Z-stack around 200 µm). All pictures were 
processed using ImageJ software (NIH, Bethesda, MD, www.nih.gov). 
Quantitative data analysis. Quantitative analysis of the whole gland reconstructions was 
performed for 10 glands from 5 mice at 8 weeks of age, 5 glands from 3 mice at 5 weeks 
of age, and 3 glands from 2 mice at 8 weeks of age (traced from 5 weeks of age). Clonal 
analysis and modelling was based on 606 subclones (157 basal subclones and 449 
luminal subclones) from 4 glands from mice at 8 weeks of age. Data was collected at 
random and all glands induced at a clonal level were included. Subclones were defined 
as the density of epithelial cells of the same type (basal or luminal) and confetti colour in 
a given branch of given level. Although our theoretical description models the distribution 
of the number of MaSCs in a TEB of level l, as we cannot access this quantity directly 
experimentally, we use as a proxy the density of labelled cells of a given type and confetti+ 
colour, in the corresponding branch of level l (i.e. the branch that was formed by the TEB 
considered in the model). Taking the density instead of the absolute number of labelled 
cells allowed us to correct for the stochastic variation of branch length that we observed 
(Extended Data Fig. 1e).  
Full Methods and any associated references are available in the supplementary 
information. 
 
  
16 
 
 
References 
 
1 Williams, J. M. & Daniel, C. W. Mammary ductal elongation: differentiation of 
myoepithelium and basal lamina during branching morphogenesis. Dev Biol 97, 
274-290 (1983). 
2 Silberstein, G. B. & Daniel, C. W. Glycosaminoglycans in the basal lamina and 
extracellular matrix of serially aged mouse mammary ducts. Mechanisms of ageing 
and development 24, 151-162 (1984). 
3 Giraddi, R. R. et al. Stem and progenitor cell division kinetics during postnatal 
mouse mammary gland development. Nature communications 6, 8487 (2015). 
4 Visvader, J. E. & Stingl, J. Mammary stem cells and the differentiation hierarchy: 
current status and perspectives. Genes & development 28, 1143-1158 (2014). 
5 Inman, J. L., Robertson, C., Mott, J. D. & Bissell, M. J. Mammary gland 
development: cell fate specification, stem cells and the microenvironment. 
Development 142, 1028-1042 (2015). 
6 Stingl, J. et al. Purification and unique properties of mammary epithelial stem cells. 
Nature 439, 993-997 (2006). 
7 Shackleton, M. et al. Generation of a functional mammary gland from a single stem 
cell. Nature 439, 84-88 (2006). 
8 dos Santos, C. O. et al. Molecular hierarchy of mammary differentiation yields 
refined markers of mammary stem cells. Proceedings of the National Academy of 
Sciences of the United States of America 110, 7123-7130 (2013). 
9 Van Keymeulen, A. et al. Distinct stem cells contribute to mammary gland 
development and maintenance. Nature 479, 189-193 (2011). 
10 Rios, A. C., Fu, N. Y., Lindeman, G. J. & Visvader, J. E. In situ identification of 
bipotent stem cells in the mammary gland. Nature 506, 322-327 (2014). 
11 Prater, M. D. et al. Mammary stem cells have myoepithelial cell properties. Nat 
Cell Biol 16, 942-950, 941-947 (2014). 
12 Bai, L. & Rohrschneider, L. R. s-SHIP promoter expression marks activated stem 
cells in developing mouse mammary tissue. Genes & development 24, 1882-1892 
(2010). 
13 Plaks, V. et al. Lgr5-expressing cells are sufficient and necessary for postnatal 
mammary gland organogenesis. Cell reports 3, 70-78 (2013). 
14 de Visser, K. E. et al. Developmental stage-specific contribution of LGR5(+) cells 
to basal and luminal epithelial lineages in the postnatal mammary gland. The 
Journal of pathology 228, 300-309 (2012). 
15 van Amerongen, R., Bowman, A. N. & Nusse, R. Developmental stage and time 
dictate the fate of Wnt/beta-catenin-responsive stem cells in the mammary gland. 
Cell stem cell 11, 387-400 (2012). 
16 Wang, D. et al. Identification of multipotent mammary stem cells by protein C 
receptor expression. Nature 517, 81-84 (2015). 
17 Zeng, Y. A. & Nusse, R. Wnt proteins are self-renewal factors for mammary stem 
cells and promote their long-term expansion in culture. Cell stem cell 6, 568-577 
(2010). 
17 
 
 
18 Badders, N. M. et al. The Wnt receptor, Lrp5, is expressed by mouse mammary 
stem cells and is required to maintain the basal lineage. PLoS One 4, e6594 
(2009). 
19 Metzger, R. J., Klein, O. D., Martin, G. R. & Krasnow, M. A. The branching 
programme of mouse lung development. Nature 453, 745-750 (2008). 
20 Sternlicht, M. D., Kouros-Mehr, H., Lu, P., and Werb, Z. (2006). Hormonal and 
local control of mammary branching morphogenesis. Differentiation, 74(7), 365- 
381. 
21 Nelson, C. M., Vanduijn, M. M., Inman, J. L., Fletcher, D. A. & Bissell, M. J. Tissue 
geometry determines sites of mammary branching morphogenesis in organotypic 
cultures. Science 314, 298-300 (2006). 
22 Lu, P. & Werb, Z. Patterning mechanisms of branched organs. Science 322, 1506-
1509 (2008). 
23 Plaks, V. et al. Adaptive Immune Regulation of Mammary Postnatal 
Organogenesis. Developmental cell 34, 493-504 (2015). 
24 Snippert, H. J. et al. Intestinal crypt homeostasis results from neutral competition 
between symmetrically dividing Lgr5 stem cells. Cell 143, 134-144 (2010). 
25 Shehata, M., van Amerongen, R., Zeeman, A. L., Giraddi, R. R. & Stingl, J. The 
influence of tamoxifen on normal mouse mammary gland homeostasis. Breast 
cancer research : BCR 16, 411 (2014). 
26  Grün, D., et al. De novo prediction of stem cell identity using single-cell 
transcriptome data. Cell Stem Cell 19, 266 – 277 (2016). 
27 Kouros-Mehr, H. and Werb, Z. Candidate regulators of mammary branching 
morphogenesis identified by genome-wide transcript analysis. Dev Dyn 12, 3404 
– 3412 (2006). 
28 Gjorevski, N. and Nelson, C.M. Integrated morphodynamic signalling of the 
mammary gland. Nat Rev Mol Cell Biol 12, 581 – 593 (2011). 
29 Ritsma, L. et al. Intestinal crypt homeostasis revealed at single-stem-cell level by 
in vivo live imaging. Nature 507, 362-365 (2014). 
30 Lopez-Garcia, C., Klein, A. M., Simons, B. D. & Winton, D. J. Intestinal stem cell 
replacement follows a pattern of neutral drift. Science 330, 822-825 (2010). 
31 Alieva, M., Ritsma, L., Giedt, R.J., Weissleder, R., & Van Rheenen, J. Imaging 
windows for long-term intravital imaging. IntraVital: 3, 2 (2014). 
32 Muraro, J., et al. A single-cell transcriptome atlas of the human pancreas. Cell 
Systems 3, 1 – 10 (2016).  
33 Hashimshony T., et al. CEL-Seq2: sensitive highly-multiplexed single-cell RNA-
Seq. Genome Biology 16, 77 (2016). 
34 Li, H. and Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler 
transform. Bioinformatics 14, 1754-1760 (2009). 
35 Anders, S. and Huber, W. Differential expression analysis for sequence count data. 
Genome Biology 11, R106 (2010). 
  
 
 
18 
 
 
Acknowledgements 
The authors would like to thank all members of the van Rheenen group for critically 
reading this manuscript, Anko de Graaff and the Hubrecht Imaging Centre for imaging 
support and the Hubrecht Institute animal caretakers for animal support. This work was 
supported by an ERC consolidator grant (648804); research grants from the Dutch 
Organization of Scientific Research (NWO; 823.02.017), the Dutch Cancer Society (KWF; 
HUBR 2009-4621), the Association for International Cancer Research (AICR; 13-0297) 
(all J.v.R), the Wellcome Trust (grant number 098357/Z/12/Z; BDS and 110326/Z/15/Z; 
EH) and equipment grants (175.010.2007.00 and 834.11.002) from the Dutch 
Organization of Scientific Research (NWO). EH is funded by a Junior Research 
Fellowship from Trinity College, Cambridge University, a Sir Henry Wellcome Fellowship 
from the Wellcome Trust and acknowledges the Bettencourt-Schueller Young 
Researcher Prize for support. CLGJS is funded by a Boehringer Ingelheim Fonds PhD 
Fellowship. 
 
Author contributions 
J.v.R, C.L.G.J.S, and A.Z. conceived the study and designed the experiments. C.L.G.J.S, 
A.Z., and N.S.M.L. performed experiments. C.L.G.J.S., E.H., M.J.M. and A.Z. performed 
analyses. C.L.G.J.S, E.H. and M.J.M made the figures. A.v.O, J.v.R. and B.D.S. 
supervised the study. All authors discussed results and participated in preparation of the 
manuscript. 
 
Competing financial interests 
The authors declare no competing financial interests. 
Corresponding authors 
Correspondence and request for materials should be addressed to J.v.R. 
(j.vanrheenen@hubrecht.eu) or B.D.S. (bds10@cam.ac.uk).  
19 
 
 
  
20 
 
 
  
21 
 
 
Figure 1 | Pubertal mammary development follows a stochastic growth pattern 
leading to heterogeneity in the ductal network. a, Confocal image of the 4th mammary 
gland showing proliferative cells (EdU+, red) exclusively in the TEBs at 3.5 weeks of age. 
Blue represents nuclei (DAPI). Dashed lines show the outline of the ducts. Scale bar 
300µm. b, Confocal image of a whole mount 4th mammary gland at 8 weeks of age stained 
for K14 (white) and DAPI (blue) alongside its manual reconstruction. Scale bar 5mm. c, 
Schematic representation of the mammary gland depicted in (b). d, TEBs either terminally 
differentiate with probability q, or bifurcate into two TEBs with probability 1-q,. e, 
Experimental probability q from reconstructions (n = 10 glands), used to predict (f). Dotted 
line indicate a balanced value of 0.5  f, Inter-subtree heterogeneity (data in purple from n 
= 10 glands) can be predicted quantitatively by a model of equipotent TEBs making 
stochastic decisions (black line, see Supplementary Theory). Left: Subtree persistence, 
defined as the distribution of subtrees having a maximal branch level i. Right: Subtree 
size distribution, where subtree size is defined as its total number of branches. Error bars 
represent mean and s.d.  
 
  
22 
 
 
  
23 
 
 
Figure 2 | Number, localization and molecular characterization of pubertal MaSCs. 
a, Reconstruction of a whole mount 5th mammary gland. b, Zoom image of the red boxed 
area in (a) showing part of a RFP+ basal clone. Scale bar 10µm. c, Schematic 
representation of the RFP+ basal clone depicted in (a) from the start of tracing until the 
most recently formed ducts. Branches containing at least one RFP+ cell are represented 
as red lines and branches without any labelled cells as black lines. d, Intravital image of 
a 5 week-old duct (left) or a TEB (right) in the 4th gland of a R26-CreERT2;R26-mTmG 
mouse . Zoom images show migration or proliferation over time (Supplementary Video). 
Scale bars 10µm. e, Top and middle: Example of the original image and outline of a TEB 
and adjacent ductal region containing a RFP+ MaSC. Bottom: An average sized 
schematic TEB (black dashed line), in which each dot represents the localization of a 
labelled MaSC in a different TEB (n = 47). f, Quantification of the number of functional 
MaSCs based on the labelling density for both the basal and luminal cell lineages. Data 
points represent different clones and symbols represent different mice (n = 3). g, t-
distributed stochastic neighbour embedding (t-SNE) map of transcriptome similarities 
between individual ductal and TEB cells isolated from a 5 weeks-old mammary gland. 
Clusters with similar transcriptional programs are highlighted with different numbers and 
colours, origin of the cells are indicated by squares (ductal cells) or ovals (TEB cells).  
  
24 
 
 
  
25 
 
 
Figure 3 | Border, but not tip MaSCs contribute short-term to ductal extension. a, 
Representative images of TEBs and the adjacent ductal region in a 5th mammary gland 
at 5 weeks of age containing confetti+ tip MaSCs (GFP, upper panel) or confetti+ border 
MaSCs (RFP, lower panel). Scale bars 100µm. b, Quantification of the clone size in the 
duct adjacent to a TEB containing either confetti+ tip MaSCs (Tip, n=10 TEBs) or confetti+ 
border MaSCs (Border, n=14 TEBs). P=0.0002, (Welch’s t-test). c, Quantification of the 
confetti clone size in the ductal region adjacent to TEBs related to the number of confetti+ 
border MaSCs (left graph) or the number of confetti+ tip MaSCs (right graph). * : P< 0.05, 
**: P< 0.01,  ***: P< 0.001 (Mann Whitney test). d, Representative image showing EdU 
labelling in a TEB and adjacent ductal region after a 3 day EdU pulse-chase experiments. 
Scale bar 100µm. e, Quantification of EdU-labelling intensity in the tip and border region 
of TEBs. P=0.03 (two-tailed t-test). Error bars represent mean and s.e.m. 
 
  
26 
 
 
  
27 
 
 
Figure 4 | Pubertal mammary gland development is driven by an equipotent pool of 
MaSCs. a, Normalized mean subclone size as a function of branch level. Subclone sizes 
were normalized by their initial value (at segment level 6). The black line represents the 
outcome of the theoretical model in the limit that cells mix completely between two 
consecutive bifurcations events. b, Clone size distributions for luminal (left) and basal 
(right) clones at successive branch levels, in log scale, each rescaled by the average 
clone size. The data displays scaling, collapsing onto an exponential dependence (black 
line) as predicted for an equipotent population. c, Whole-gland reconstruction, with the 
fraction of confetti colours present in branches represented as pie charts for each level.  
d, Representative confocal images of ducts and TEBs from different areas of the 4th 
mammary gland depicted in Extended Data Fig 11. RFP+ cells in red, GFP+ cells in green 
and DAPI in blue. Scale bars 100µm. Error bars represent mean and s.e.m. 
  
28 
 
 
  
29 
 
 
Figure 5 | Model of pubertal mammary ductal growth dynamics. Cartoon depicting 
the proposed model of pubertal mammary ductal outgrowth: Highly proliferative MaSCs 
are located in the TEB. Progeny of the border MaSCs is left behind as the TEB is pushed 
forward due to proliferation. Early ductal cells lose stemness and undergo a defined 
number of cell divisions to produce terminally differentiated ductal cells (cluster). During 
the elongation process, MaSCs randomly mix in the TEB. At the start of the bifurcation 
process, MaSCs are doubled and randomly partitioned between daughter TEBs leading 
to clonal enrichment and extinction. Subsequently, a fraction of TEB terminate, with their 
MaSCs all stopping contributing to growth, while another fraction repeats another cycle 
of branching and elongation. 
  
30 
 
 
Extended Data 
“Identification of mammary stem cells and their dynamics during branching 
morphogenesis” by Scheele et al. 
 
Extended Data Figure 1 | Pubertal mammary gland development is driven by elongation  
and bifurcation of proliferative TEBs.          32 
 
Extended Data Figure 2 | Unbiased lineage tracing strategy during pubertal development.   36 
 
Extended Data Figure 3 | Comparisons between branch patterns at different doses  
of tamoxifen and of 4th vs 5th mammary glands.       38 
 
Extended Data Figure 4 | Confetti+ mammary glands reveal the contribution of MaSCs during  
pubertal development.          40 
 
Extended Data Figure 5 | MaSCs are localized in the TEBs and contribute during oestrous driven  
proliferation and alveoli formation during pregnancy.      42 
 
Extended Data Figure 6 | Clonal labelling enables quantification of the number of functional  
MaSCs based on labelling density.         44 
 
Extended Data Figure 7 | StemID identifies basal and luminal cell clusters in a mixture of  
mammary epithelial cells.          46 
 
Extended Data Figure 8 | Successive rounds of bifurcation lead to monoclonal conversion  
as pubertal development proceeds, at a speed dictated by the rate of MaSC mixing in TEBs.   48 
 
Extended Data Figure 9 | Successive rounds of bifurcation lead to clonal extinction and  
31 
 
 
enrichment.           50  
 
Extended Data Figure 10 | Theory of equipotent TEBs making near-balanced stochastic fate  
choices between branching and termination explains quantitatively kidney branching  
morphogenesis.           52 
  
32 
 
 
  
33 
 
 
Extended Data Figure 1 | Pubertal mammary gland development is driven by 
elongation and bifurcation of proliferative TEBs. a, Confocal image of a whole-mount 
mammary glands at the onset of puberty (3.5 weeks of age) stained for DAPI (in blue) 
and EdU (magenta, 4 hours chase). Ducts are outlined in yellow. Scale bar represents 5 
mm. b, Zoom of part of the 5th mammary gland outlined with a red box in (a). Confocal 
image shows a single Z-plane of a proliferative TEB in the 5th gland at 3.5 weeks of age. 
Scale bar represents 50 µm. c, Left: Zoom of the 5th mammary gland outlined with a 
yellow box in (a). Top panel shows the maximum projection of the 5th gland at 3.5 weeks 
of age. Bottom panel shows a zoom of a representative TEB (single Z-plane, outlined with 
a white box). Middle: Maximum projection of a 5th mammary gland at 5 weeks of age 
Bottom panel shows a zoom of a representative TEB (outlined with a white box). Right: 
Maximum projection of the 5th mammary gland at 8 weeks of age. Bottom panel shows a 
zoom of a representative TEB (single Z-plane, outlined with a white box). EdU positive 
cells shown in magenta, nuclei are stained with DAPI (in blue), yellow outlines mark the 
ducts. Scale bars represent 300 µm. d, Schematic representation of the methods used to 
count the branch levels in the mammary gland. In pubertal developing mammary glands, 
younger ducts (light grey) that were formed more recently have undergone more 
bifurcation events than older ducts (dark grey) that were formed earlier. e, Cumulative 
Distribution of branch length at 5 weeks and 8 weeks, for branches formed after the 
beginning of puberty (level > 6). f, Bar graph showing the fraction of EdU+ cells in 
proliferative end buds during the course of pubertal development (n =3 mice per time 
point, 10 TEBs per mouse). g, At 3.5, 5, and 8 weeks (n = 3 mice per condition), the 
fraction of EdU+ TEBs were measured. The graph shows the fraction of EdU+ TEB plotted 
against the highest branch level. The black line represents the prediction from the 
stochastic model of equipotent TEBs described in the main text. h, Probability q for a TEB 
to terminally differentiate as a function of branch level, read from the whole-gland 
reconstructions (n = 10 glands, same data as Fig. 1e), overlaid with an linear piecewise 
fit (dashed green lines) which we use as an input to predict the subtree heterogeneity in 
Fig. 1f. Dotted lines indicate respectively a constant value of 0.5 (horizontal line), 
indicative of balance between TEB bifurcation and termination, and the branch level at 
34 
 
 
the start of puberty (vertical line). i, Normalized average branch width as a function of 
branch level. Dotted lines indicate resp. a constant value of 1 (horizontal), and the branch 
level at the start of puberty (vertical). j, Normalized average branch length as a function 
of branch level. Dotted lines indicate respectively a constant value of 1 (horizontal), and 
the branch level at the start of puberty (vertical). k, Fate correlation between related TEBs, 
as a function of branch level. Correlation between sibling TEBs is in purple, while 
correlation between first-cousin TEBs is in green. A value of 1 indicates that the fate of 
two related TEBs is perfectly correlated, whereas a value of -1 indicates perfect anti-
correlation. The dotted line indicates 0, i.e. no correlation. Error bars represent mean and 
s.e.m. For all graphs in this figure, data of the 4th and 5th mammary gland were combined. 
 
  
35 
 
 
  
36 
 
 
  
37 
 
 
Extended Data Figure 2 | Unbiased lineage tracing strategy during pubertal 
development. a, Cartoon depicting the gene system used for the lineage tracing strategy. 
b, Timeline of the lineage tracing strategy. c, Labelled cell fraction (product of clonal 
persistence and the surviving mean subclone size) as a function of segment level, for the 
3 to 8 weeks tracing. The labelled cell fraction hovers around its initial value (reference 
taken for segment 6), indicative of a representative labelling. This also signifies that the 
persistence decays inversely to the linear increase in clone size, as expected from the 
model of long-term equipotent MaSCs. d, Labelled cell fraction as a function of segment 
level, for the 5 to 8 weeks tracing. e, Labelled cell fraction as a function of segment level, 
for the 3 to 5 weeks tracing. In both cases, the labelled cell fraction hovers around its 
initial value, confirming the results of the 3 weeks to 8 weeks of age tracing. f, Normalized 
mean subclone size (basal in purple and luminal in green) as a function of branch level, 
for the 3 to 5 weeks tracing, showing a linear increase as branching morphogenesis 
proceeds due to cell segregation upon bifurcation. Subclone size is normalized by the 
area of the segment they are located on, in order to correct for the segment length 
variability. Subclone size was also normalized to its initial value (at segment level 6), and 
display a similar trend as the 3 to 8 weeks tracing. g, Clone size distributions for luminal 
(left) and basal (right) clones at 5 weeks of age (lineage tracing from 3 weeks to 5 weeks 
of age), at successive branch levels, in log plot. The clone size was rescaled by the 
average clone size at each time point. In each case, the data displays scaling, and 
collapses on a single-exponential (black line), as expected for an equipotent cell 
population. h, Quantification of the number of basal and luminal cells per µm in mammary 
ducts of the adult mammary gland (n = 17 randomly picked ducts from three different 
mice). Error bars represent mean ± s.e.m. For all graphs in this figure, data of the 4th and 
5th mammary gland were combined. 
  
38 
 
 
39 
 
 
Extended Data Figure 3 | Comparisons between branch patterns at different doses 
of Tamoxifen and of 4th vs 5th mammary glands. a, Confocal images of whole-mount 
4th mammary glands at 8 weeks of age, after injection of 3mg/25g, 1.5mg/25g, or no 
Tamoxifen respectively at the onset of puberty (3 weeks of age). Injection of 3mg/25g 
clearly interferes with branching morphogenesis, as the mammary gland did not fully 
invade the fat pad at the end of puberty. The mammary gland treated with 1.5mg/25g 
does not show any delay in the mammary gland development, and shows the same 
branching pattern as the control gland. b, Normalized average branch width and length 
as a function of branch level. No differences were observed in branch length or branch 
width for the different concentrations of Tamoxifen. c, Probability q for a TEB to terminally 
differentiate as a function of branch level, read from the whole-gland reconstructions. In 
b and c, we used resp. n = 2 glands (no Tamoxifen), n = 2 glands (0.2mg/25g), and n = 
5 glands (1.5mg/25g). No differences were observed between Tamoxifen induced glands 
and the control gland. For graphs b,c, data of the 4th and 5th mammary gland were 
combined. d, Normalized average branch widths and lengths as a function of branch 
level. e, Probabilities q for a TEB to terminally differentiate as a function of branch level, 
read from the whole-gland reconstructions. f, Cumulative subtree size distributions, where 
size is defined as the total number of branches in a given subtree (starting at level 6). g, 
Subtree persistence, defined as the distribution of subtrees (starting at level 6) having a 
given maximal branch level i. In all cases, we considered 4th glands (n = 5) and 5th glands 
(n = 5), and plot comparisons respectively in purple and in green. Error bars represent 
mean and s.e.m (b-e) and s.d. (f,g). 
  
40 
 
 
41 
 
 
Extended Data Figure 4 | Confetti+ mammary glands reveal the contribution of 
MaSCs during pubertal development. a, Confocal image of a whole-mount confetti+ 4th 
mammary gland at 5 weeks of age, traced from 3 weeks of age. The ducts were manually 
reconstructed by outlining the ducts based on the K14 staining (white) and all individual 
confetti-labelled cells were outlined with their respective confetti colour (blue symbols: 
CFP+ cells; green symbols: GFP+ cells; yellow symbols: YFP+ cells; red symbols: RPF+ 
cells). b, Image shows a duct in the rudimentary gland containing single confetti labelled 
cells (outlined with the red box in (a)). c, d, Confocal images showing ducts that were 
formed after the induction of confetti-tracing containing clonal outgrowths. Notice that only 
one confetti colour is present in these ducts close to the border of the fat pad. Scale bars 
represent 100 µm. e, Probability for a clone to be lost from the subsequent subtree, as a 
function of its last segment level. f, Distribution of intermittency length (purple dots), 
defined as the number of consecutive unlabelled branches existing between branches 
labelled in the same colour and cell type. The black line represents an exponential fit. 
Error bars represent mean ± s.e.m. 
  
42 
 
 
 
43 
 
 
Extended Data Figure 5 | MaSCs are localized in the TEBs and contribute during 
oestrous driven proliferation and alveoli formation during pregnancy. a, Cartoon 
depicting the gene system used for the lineage tracing strategy and timeline of the lineage 
tracing strategy. b, Confocal images of a duct and a lobuloalveolar structure in the 4th 
mammary gland containing lineage traced cells after 18 months of tracing contributing to 
ductal turnover and the formation of lobuloalveolar structures driven by the oestrous 
cycle. c, Cartoon depicting the gene system used for the lineage tracing strategy and 
timeline of the lineage tracing strategy. d, Maximum projection of a part of a whole mount 
4th mammary gland after 4 days of lactation, showing the contribution of YFP+ confetti 
cells to alveoli formation. Scale bars represent 100 µm. e, Whole mount image of a 5th 
mammary gland at 5.5 weeks of age, traced from 3 weeks of age. Ducts are stained for 
K14 (in red). Scale bar represents 5mm. f, Maximum projection of a part of the 5th 
mammary gland depicted in (e) showing a luminal GFP clone contributing to the subtree 
depicted here. Scale bar represents 1mm. g-h, Maximum projection of TEBs and their 
adjacent ducts outlined by the yellow and green boxes in (f). The most distal GFP cells of 
the clone are situated in the TEBs, and are actively dividing as indicated by the most 
recent progeny located in the adjacent ducts. The localization of these most distal cells 
of a clone was determined as shown in Fig. 2e. Scale bar represents 100 µm. 
  
44 
 
 
  
45 
 
 
Extended Data Figure 6 | Clonal labelling enables quantification of the number of 
functional MaSCs based on labelling density. a, Cartoon depicting the method used 
to determine the number of MaScs. Clonal labelling at the onset of puberty results in the 
confetti labelling of 1 MaSC out of N MaSCs in a TEB (red cell). Clonal analysis in the 
resulting subtree at the end of puberty enables to determine the contribution of 1 labelled 
MaSC to the entire subtree that is formed after the induction of confetti-tracing. The 
labelled fraction is calculated by dividing the clone size (red cells) through the total 
number of cells in the subtree (1/N). By taking the inverse, the number of functional 
MaSCs N can be calculated. b, Representative images of a part of a 4th whole mount 
mammary gland, stained for K14 (basal marker) used to determine the identity of the 
labelled clones. Zooms show a single Z-plane with basal cells (red outline) which overlap 
with the staining, or luminal cells (yellow outline) which do not overlap with the basal 
staining. Scale bars represent 100 µm. c, Maximum projection of part of a 5th whole mount 
mammary gland, stained for E-cadherin (luminal marker) used to determine the identity 
of the labelled clones. Zooms show a single Z-plane with basal cells (yellow outline), 
which do not overlap with the staining, or luminal cells (red outline), which do overlap with 
the luminal staining. Scale bar left panel is 1 mm, scale bars right panels is 100 µm. d, 
Quantification of number of basal and luminal cells in the TEB (n = 10 TEBS from 3 
different mice). Bar graphs represent mean ± s.e.m. e, Representative images of TEBs 
at 5 weeks of age showing a single Z-plane of the TEB stained for DAPI and K14 or E-
cadherin to determine the number of luminal and basal cells per TEB. Z-stacks of 
complete TEBs (Z-step size 5 µm) were used to count the total number of basal cells per 
TEB. Scale bars represent 100 µm.  
  
46 
 
 
  
47 
 
 
Extended Data Figure 7 | StemID identifies basal and luminal cell clusters in a 
mixture of mammary epithelial cells. a, Heatmap of cell-to-cell transcriptome distances 
measured by 1- Pearson’s correlation coefficient. Basal and luminal cell clusters identified 
by StemID are colour coded along the axes. b, Cumulative transcript counts (colour 
legend) of cell cycle genes expressed during G1 phase, S phase and M phase 
respectively are indicated in the t-SNE maps. None of the cell cycle gene groups show a 
clear correlation with any of the clusters. c, Transcript counts (colour legend) of basal cell 
markers smooth muscle actin, K5 and K14 are indicated in the t-SNE maps. Expression 
of these genes is restricted to clusters 6, 7, 8, and 9. d, Transcript counts (colour legend) 
of luminal cell markers K8, K18 and K19 are indicated in the t-SNE maps. Expression of 
these genes is restricted to clusters 1, 2, 3, and 5. e, Transcript counts (colour legend) of 
previously reported MaSC markers Axin2 and Lgr5 are indicated in the t-SNE maps. Left 
graph: Expression of Axin2 is assigned to the basal cells, both ductal and TEB cells, but 
is not restricted to one identified StemID cluster. Right graph: Expression of Lgr5 is 
assigned to both luminal and basal cell clusters, and is not restricted to one identified 
StemID cluster. Note that the transcript counts for both genes are very low. f, g, Heatmap 
of differentially expressed genes comparing luminal cell cluster 1 with luminal cell clusters 
2 and 3 (f), and comparing basal cell clusters 7 and 8 with basal cell clusters 6 and 9 (g). 
Rows are genes and columns are cells ordered based on cluster number. Genes are 
grouped based on hierarchical clustering. Log10 expression of transcript counts is plotted.  
  
48 
 
 
49 
 
 
Extended Data Figure 8 | Successive rounds of bifurcation lead to monoclonal 
conversion as pubertal development proceeds, at a speed dictated by the rate of 
MaSC mixing in TEBs.  a, Simulated lineage tree in two opposite limits. Top: Well-mixed 
limit, where cells are randomly segregated upon bifurcation, showing slow linear 
conversion towards monoclonality. Bottom: Zero-mixing limit, where the position of cells 
during two successive bifurcations remains unchanged, showing a fast exponential 
conversion towards monoclonality. Branches without any labelled cells are represented 
as black lines. Branches containing confetti-labelled cells are represented as lines of the 
dominant colour, and branching points show a pie chart representing the proportion of 
each colour present in the branch. The size of the pie chart represents the size of the 
confetti clone to the total duct. b, Mean persisting confetti-clone size as a function of 
branch level, in the same two opposite limits, assuming a total MaSC number of N = 200. 
Numerical integration in the well-mixed limit (top), for three initial conditions for the initial 
number of cells labelled (n = 1 in blue, n = 2 in yellow, n = 3 in green). The red dashed 
line indicates the predicted linear slope. Numerical integration in the zero-mixing limit 
(bottom, n = 1 in blue). The red dashed line indicates the predicted exponential increase 
for small clone sizes. c, Top: Confocal images of TEBs during the bifurcation process. 
Bottom: Schematic representation of the branching process as predicted by the model in 
the well-mixed limit. d, Evolution of the cumulative rescaled subclone size distribution at 
different segment numbers (7,9,11,15 resp. in purple, green, blue and orange) for a two-
compartment model of TEB, consisting of a 10% population of SCs tilted towards renewal 
and 90% population of progenitors tilted towards loss. Contrary to the equipotent one-
population model and our data of Fig. 3, the two-compartment model does not display 
scaling, with SCs forming a long-tail in the distribution. e, Representative confocal images 
of a confetti+ mammary epithelial duct of a 4th mammary gland at 8 weeks of age. Clones 
of RFP+ cells (red) are outlined and blue represents E-cadherin. TEBs are located at the 
right side of the image. Scale bars represent 300µm. f, Quantification of the average 
cluster size for basal and luminal cell lineages at 8 weeks of age. P=0.46 (Mann-Whitney 
test). 
  
50 
 
 
51 
 
 
Extended Data Figure 9 | Successive rounds of bifurcation lead to clonal extinction 
and enrichment. a, Confocal image of a whole mount-imaged 4th mammary gland 
stained for K14 (white) and DAPI (blue). The ductal network was manually reconstructed 
by outlining the ducts based on the K14 staining (white) and all individual confetti-labelled 
cells were outlined with their respective confetti colour (blue symbols: CFP+ cells; green 
symbols: GFP+ cells; yellow symbols: YFP+ cells; red symbols: RPF+ cells). This 
information was used as input for a schematic representation of the lineage tree in Fig. 
3d (lower panel). Scale bar represents 5 mm. b-d, Confocal images showing a zoom of 
the location of the colored boxes indicated in (a), showing a conversion of these luminal 
clones towards monoclonality over successive bifurcation rounds, including extinction 
and enrichment of clones. Scale bars represent 100 µm. 
  
52 
 
 
 
53 
 
 
Extended Data Figure 10 | Theory of equipotent TEBs making near-balanced 
stochastic fate choices between branching and termination explains quantitatively 
kidney branching morphogenesis. a, Maximum projection of a whole mount image of 
a murine kidney at E16. Ducts are stained for E-Cadherin (green) and nuclei for DAPI 
(blue). Scale bar represents 1mm. b, Outline from our quantitative reconstruction. c, 
Probability q for a bud to terminally differentiate as a function of branch level, read from 
the whole-gland reconstruction. Green line represents a constant value of 0.5, indicative 
of balance between TEB bifurcation and termination. d, Distribution of branch lengths 
(level > 6), which collapses onto an exponential distribution, indicative of branching 
occurring stochastically at any time with equal constant probability. e,f, Inter-subtree 
heterogeneity (level >6, data in purple) can be predicted quantitatively by a model of 
equipotent buds making stochastic decisions: subtree persistence (e), defined as the 
distribution of subtrees having a given maximal branch level i, and subtree size 
distribution (f), where size is defined as the total number of branches in a given subtree. 
The theoretical curve (green line) explains well the subtree heterogeneity on both 
measurements. Error bars represent mean and s.e.m (c) or mean and s.d. (d,e). 
 
 
  
54 
 
 
Supplementary Information 
“Identification of mammary stem cells and their dynamics during branching 
morphogenesis” by Scheele et al. 
 
Supplementary Methods        55  
     
  
55 
 
 
Supplementary Methods 
 
Whole mount immunofluorescence staining of mammary glands. The fourth and the 
fifth mammary glands were dissected and incubated in a mixture of collagenase I (1mg/ml, 
Roche Diagnostics) and hyaluronidase (50 μg/ml, Sigma Aldrich) at 37°C for optical 
clearance, fixed in periodate-lysine-paraformaldehyde (PLP) buffer (1% 
paraformaldehyde (PFA, Electron Microscopy Science), 0.01 M sodium periodate, 0.075 
M L-lysine and 0.0375 M P-buffer (0.081 M Na2HPO4 and 0.019 M NaH2PO4) (pH 7.4)) 
for 2 hours at room temperature (RT), and incubated for 2 hours in blocking buffer 
containing 1% bovine serum albumin (Roche Diagnostics), 5% normal goat serum 
(Monosan) and 0.8% Triton X-100 (Sigma-Aldrich) in PBS. Primary antibodies were 
diluted in blocking buffer and incubated overnight at RT. Secondary antibodies diluted in 
blocking buffer were incubated for at least 4 hours. Nuclei were stained with 4',6-
diamidino-2-phenylindole (DAPI, 0.1 μg/ml; Sigma-Aldrich) in PBS. Glands were washed 
with PBS and mounted on a microscopy slide with Vectashield hard set (H-1400, Vector 
Laboratories). Primary antibodies: anti-K14 (rabbit, Covance, PRB155P, 1:700) or anti-
E-cadherin (rat, eBioscience, 14-3249-82, 1:700). Secondary antibodies: goat anti-rabbit 
or goat anti-rat, both conjugated to Alexa-647 (Life Technologies, A21245 and A21247 
respectively, 1:400). For EdU cell proliferation staining of whole mount mammary glands, 
a click-it staining (Click-iT™ EdU, Invitrogen) was performed according to the 
manufacturer’s instructions before staining with the primary antibodies as described 
above. 
56 
 
 
Whole mount immunofluorescence staining of embryonic kidney. Kidneys were 
dissected from embryos at day E15.5, and fixed in PLP buffer overnight at 4°C, and 
incubated for 4 hours in blocking buffer containing 1% bovine serum albumin (Roche 
Diagnostics), 5% normal goat serum (Monosan) and 0.8% Triton X-100 (Sigma-Aldrich) 
in PBS. Primary antibodies were diluted in blocking buffer and incubated overnight at RT. 
Secondary antibodies diluted in blocking buffer were incubated for at least 6 hours. Nuclei 
were stained with 4',6-diamidino-2-phenylindole (DAPI, 0.1 μg/ml; Sigma-Aldrich) in PBS. 
 
EdU incorporation experiments. 3.5, 5 or 8 week-old mice (n = 3 mice for each age) 
were injected IP with 0.5 mg 5-ethynyl-2′-deoxyuridine (EdU, Invitrogen) diluted in 
phosphate buffered saline (PBS). For the EdU pulse-chase experiments, EdU (0.5 mg) 
was IP injected in 4 week-old mice (n = 3). Mice were sacrificed 4 or 72 hours after EdU 
injection and the fourth and fifth mammary glands were collected and processed as whole 
mount glands. For analysis, three-dimensional tile scan images of the whole mount 
glands were taken and the number of EdU+ TEBs was scored. For the quantification of 
the fraction of EdU+ cells, 10 EdU+  TEBs per time point were selected in a blinded manner 
and the number of EdU+ cells were counted manually for each selected TEB. For the 
pulse-chase experiment, 3 EdU+ TEBs per mouse were selected in a blinded manner and 
the intensity of 10 randomly picked EdU+ cells was measured for both the tip and the 
border area of the selected TEB and a two-tailed t-test was performed. Normal distribution 
was confirmed using a d’Agostino and Pearson omnibus normality test. The variance 
57 
 
 
between the groups was tested with an F-test and was found to be not significantly 
different. 
 
Analysis of whole mount imaged mammary glands and embryonic kidney. Three-
dimensional tile scan images of whole mount mammary glands and embryonic kidneys 
were used to manually reconstruct the ductal network by outlining the ducts. Labelled 
confetti cells were annotated in the schematic outline of the mammary tree including 
information on the confetti colour for the mammary glands (GFP = green, YFP = yellow, 
RPF = red, and CFP = cyan). In these schematic outline images, the length and width of 
all the ducts, the coordinates of the branch points, and the position of the labelled cells in 
ducts and in TEBs were scored, which was used as input for a schematic representation 
of the lineage tree. To depict the topology of the resulting tree, the Newick format was 
used, a widely used format to represent hierarchical structures using nested parentheses 
to encode information about the linkage between branches, their respective length and 
number of cells of each confetti colour. Custom-made Python software and the ETE2 
python toolkit were used for the conversion and visualization of the schematic mammary 
gland lineage tree. The origin of the gland was always located at the top of the 
reconstruction.  
 
Intravital imaging of the mammary glands. R26-CreERT2;R26-mTmG mice were IP 
injected with 0.2mg/25gTamoxifen diluted in sunflower oil (Sigma) at 3 weeks of age. At 
58 
 
 
5 weeks of age, a mammary window was inserted near the 4th and 5th mammary gland, 
for details see ref 31. Mice were anesthetized using isoflurane (1.5% isoflurane/medical 
air mixture) and placed in a facemask with a custom designed imaging box. Mice were IP 
injected with AZD-7762 (0.5 mg in PBS, Sigma) every 5 hours during the time lapse 
imaging. Imaging was performed on an inverted Leica SP8 multiphoton microscope with 
a chameleon Vision-S (Coherent Inc., Santa Clare, CA, www.coherent.com). The 
microscope is equipped with four HyD detectors: HyD1 (<455 nm), HyD2 (455–490 nm), 
HyD3 (500–550 nm), and HyD4 (560–650 nm). Different wavelengths between 700 nm 
and 1,150 nm were used for excitation; collagen (second harmonic generation) was 
excited with a wavelength of 860 nm, and detected in HyD1. GFP and Tomato were 
excited with a wavelength of 960 nm, and detected in HyD3 and HyD4. TEBs and ducts 
were imaged every 20-30 minutes, using a Z-step size of 3 µm over a minimum period of 
8 hours. All images were in 12 bit and acquired with a 25× (HCX IRAPO N.A. 0.95 WD 
2.5 mm) water objective. 
 
FACS sorting. Single TEBs and ducts were isolated from 5 week-old MMTV-Cre;R26-
loxP-stop-loxP-YFP mice. YFP expression was used to determine the localization and 
structure of the mammary gland. Single TEBs and pieces of ducts were micro dissected 
from the 4th and 5th mammary glands using micro dissection tools. Single TEBs and ducts 
were digested in DMEM/F12 (GIBCO, Invitrogen Life Technologies) supplemented with 
hyaluronidase (300 µg/ml, Sigma Aldrich), and collagenase I (2mg/ml,Roche Diagnostics) 
at 37⁰C, followed by centrifugation at 550g for 10 minutes. The fatty layer on top and the 
59 
 
 
aqueous layer in the middle were aspirated, and the remaining pellet was dissolved in 
5mM EDTA/PBS with 5% fetal bovine serum (Sigma) and kept on ice for 10 minutes 
before labelling with the following antibodies: anti-mouse CD45-Pacific blue (clone 30-
F11, Biolegend), anti-mouse CD31-Pacific blue (clone 390, Biolegend), and rat anti-
mouse CD140a-Pacific blue (Clone APA5, BD Bioscience). Cells were incubated for 30 
minutes on ice in the dark, washed once in 5mM EDTA/PBS with 5% fetal bovine serum 
(Sigma), and centrifuged at 250g for 5 minutes. Pellet was dissolved in 5mM EDTA/PBS 
with 7-AAD, and sorted on a FACS AriaII Special Ordered Reseach Product (BD 
Biosciences). A broad FSC SSC gate was followed by a gate excluding doublets. Next, 
7-AAD negative cells were gated, and from this population Lin- (CD45-, CD31-, CD140a-) 
and YFP+ cells were sorted into 384-well plates containing 5 µl of mineral oil that 
contained a droplet of 200 nl primers, dNTPs and synthetic mRNA molecules (ERCC). 
Single cell mRNA sequencing. Single-cell mRNA sequencing was done as described 
in ref 32. Briefly, cells were sorted into 5 µl of mineral oil that contained a droplet of 200 
nl primers, dNTPs and synthetic mRNA molecules (ERCC). Cells were fused with this 
droplet by centrifugation and lysed at 65˚C, followed by RT and second strand synthesis 
aided by a Nanodrop II liquid handling platform. The resulting cDNA was processed 
following the CEL-Seq2 protocol33. Libraries were sequenced on an Illumina NextSeq 
with 75bp paired-end reads. The 5’mate was used to identify cells and libraries while the 
3’mate was aligned to the mm10 RefSeq mouse transcriptome using BWA34. Analysis 
was done using StemID, for details of the methodology see reference 26. Endothelial 
cells, erythrocytes and lymphocytes were removed from the dataset by filtering out cells 
60 
 
 
based on expression of Cd36 (5 unique transcripts), Beta-s (1000 unique transcripts) and 
Cd74 (2 unique transcripts)/Cd52 (1 unique transcript) respectively (27 cells in total). The 
remaining cells were normalized by down sampling to 3000 transcripts, after which 
StemID26 was used for clustering and cell type annotation. Cells with less than 3000 
unique transcripts were discarded. In total 91 cells were included, of which we could 
assign 36 cells to the luminal lineage (cluster 1 containing 9 cells, cluster 2 containing 17 
cells, cluster 3 containing 2 cells, and cluster 5 containing 10 cells), 51 cells to the basal 
lineage (cluster 6 containing 14 cells, cluster 7 containing 10 cells, cluster 8 containing 8 
cells, and cluster 9 containing 17 cells), and 4 cells to a non-epithelial origin (cluster 4). 
After identifying luminal and basal cell clusters, three cells that were erroneously 
annotated as basal cells belonging to cluster 8 were manually annotated to belong to 
luminal cluster 2, based on luminal/basal markers like Krt8/Krt18 and Krt5/Acta2 
respectively. Cluster 5 expressed higher levels of the pre-ribosomal RNA 45s, which is 
often found in cells of low quality. We therefore chose to omit cluster 5 from further 
analysis. Differential gene expression between clusters was based on a previous 
method35 and done as described in ref 32. All data analysis with StemID and custom 
scripts was performed with Rstudio, version 0.99.491. The data discussed in this 
publication have been deposited in NCBI’s Gene Expression Omnibus, and are 
accessible through GEO Series accession number GSE85875 
(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE85875). 
 
 
61 
 
 
 
  
62 
 
 
Supplementary Theory 
“Identification of mammary stem cells and their dynamics during branching 
morphogenesis” by Scheele et al. 
 
     
Supplementary Theory        63 
In this Supplementary Theory Note, we provide further details on our modelling
approach of pubertal mammary gland development in the mouse during puberty, as
well as on the data analysis and statistics used. We start by describing the meso-
scopic branching choices leading to the topology of the tree, before turning to the
clonal data. Throughout the note, we refer to three length scales, which we now
define: the macroscopic scale refers to the structure of the gland as a whole, the
mesoscopic scale refers to terminal end buds (TEBs) and their proximate ducts, and
the microscopic scale refers to basal and luminal cells.
1 Fate choices of terminal end buds leading to glan-
dular structure
Given the very large amount of data provided by whole-gland reconstruction coupled
to lineage tracing, the first question we need to address is whether different spatial
regions of the gland have the same statistical characteristics, or if some are predes-
tined to grow while some are predestined to differentiate early. This is crucial as it
conditions the subsequent interpretation of the lineage tracing data, by determining
at which length scale to look at average properties. In previous studies, it has in-
deed been noted and observed that the mammary gland displays a non-stereotypic
branching [36, 37], in opposition to the lung [38], which could in principle complexify
the analysis of the lineage tracing data. Nevertheless, as we state in the main text,
we show in the following sections that the non-stereotypic pattern of the mammary
gland betrays a very simple rule underlying its growth.
From the whole-gland reconstructions, we examine successively two different prop-
erties of the ductal tree organisation: its topology, i.e. how different branches are
connected to each other, without considering the geometrical properties of individual
branches, and its detailed spatial structure, i.e. how the width, length and geometry
of individual branches change in space and time throughout development.
1.1 Topology of the mammary gland ductal epithelium
Since there are no reconnections, the topology of the tree is defined unambiguously
by the behaviour of TEBs alone, which only have two choices (Fig. 1d). A TEB i may
branch into two new TEBs (with probability 1−qi), or terminate through differentiation
and stop contributing to growth (with probability qi). Therefore, the full topology of a
tree containing n segments can be deduced from n choices made by TEBs, which we
note wi(t)i=1..n, where wi = 0 corresponds to termination and wi = 1 corresponds to
branching. Moreover, since we have at our disposal whole-tree reconstructions, we
can simply read off from the final structure of the tree the ”historical” choices made
by each TEB. In principle, each of these events could be genetically programmed by
a specific signalling cascade, or could be correlated spatio-temporally among each
other. However, and as detailed in the main text, we explore first the simplest hy-
pothesis of TEBs each behaving according to the same set of probability (which may
evolve in time), and making stochastic decisions between branching and termination.
We test in the following whether this hypothesis can explain quantitatively the data,
and check its assumptions and implications experimentally. It should be noted that
there are also rare events (≈ 5%) involving TEB trifurcations, although from recon-
structions it is difficult to determine if these are “real” features, or rather two close
consecutive branching events. For the sake of consistency, we assume in the fol-
lowing that 3-way branching consisted of two consecutive branching events, with the
order taken at random.
Examining the final topology of the gland does not per se allow us to access the
exact timing at which a decision wi was taken. However, we will use in the following
the generation index, or branch level, as a proxy of time: the root of the mammary
tree, in contact with the nipple, is at level 0 (Extended Data Fig. 1d). Again, the
relevant of this proxy will be verified later in the study. The level number of a branch
segment is the minimal number of branch points one must pass to reach the tree root.
We therefore average fate choices wi over all branches of the same level number l,
and examine the evolution of this quantity as a function of segment level l.
Therefore, denoting B an active/proliferative TEB, and D a differentiated TEB, the
tree structure in the equipotent TEB hypothesis is entirely determined by the following
dynamics:
B→
{
D Proba. q(l)
B+B Proba. 1−q(l) (1)
which translates to a stochastic Galton-Watson type birth-death process.
At a given level l, or developmental time, if q(l)< 1/2, the number of active TEBs
increases, whereas for q(l)> 1/2, this number decreases. It should be noted though
that even for q(l)= 1/2 (perfect balance), the organ as a whole would continue to grow.
In this case, the number of active TEBs is conserved on average, while the number
of differentiated TEBs increases linearly with each branching event. However, most
branched organs would stochastically have all their TEBs terminating over time, while
an ever smaller subset of trees would continue growing.
Therefore, according to this dynamics, the total average number of branches S(l)
at level l is obtained by the kinetic equation, in the continuous limit,
dS
dl
= (1−2q(l))S(l) (2)
so that, starting from a single TEB at level 0, S(l) = e∑
l
0(1−2q(l))dl.
1.2 Strategies of branching morphogenesis
In principle, there are many ways of building a branched organ. TEBs could for in-
stance always branch throughout development, before differentiating collectively at
the end of the process, at a level l f (bang-bang control). Within the framework of the
dynamics above, this corresponds to the condition,
q(l) =
{
0 if l < l f
1 if l > l f
(3)
and gives rise to a deterministic, exponential growth process. Two glands from two
different animals are then expected to be identical from a topological perspective.
At the other extreme, another strategy to build a branched organ is that TEBs
branch asymmetrically at the population level, giving rise on average to an active
TEB and an inactive one, until level l f . This corresponds to
q(l) =
{
1/2 if l < l f
1 if l > l f
(4)
Within this strategy, there are again two alternative possibilities. Either TEBs have
an engrained capacity to branch asymmetrically, and always give rise to only one pro-
liferative offspring, or TEBs make the choice between branching and termination is a
purely stochastic manner, uncorrelated to the outcomes of its ancestors/neighbours.
In this case, a population level balance has to be implemented, and could be im-
plemented either intrinsically (for instance by bistable switches and negative interac-
tions between cells being segregated, leading to on average one TEB getting termi-
nated, although this would be very complex) or extrinsically (for instance because of
crowding-effects, discussed in more depth below, or because of hormonal or stromal
external signals [37]).
Interestingly, this bears very strong parallels, at the organ scale, with the two dif-
ferent strategies of fate choices that individual stem cells can make. Historically, the
consensus was of deterministic ”immortal” stem cells with engrained asymmetric di-
vision (analogous to the first hypothesis above). However, we and other groups have
shown in the past ten years that many, if not most, epithelial stem cells and progeni-
tors make stochastic fate choices at the individual level [39, 40], and that the balance
in fate is only encoded at the population level (analogous to the second hypothesis
above). In the context of individual stem cell fate, the same question of intrinsic vs
extrinsic regulation can be raised, and remains largely un-solved: are stochastic fate
probabilities encoded intrinsically within individual cells, or are they the manifestation
of differential access to extrinsic signals, such as mechanical forces or niche factors?
1.3 Fitting strategy and results of the modelling
1.3.1 Fitting the evolution of the branching probability
Turning to the dataset, we then used the whole gland reconstructions to measure
q(l), the average probability for a TEB at level l to branch (Extended Data Fig. 1h).
Strikingly, we found that:
• towards the base of the tree (l < 6), i.e. the part that pre-existed at the beginning
of puberty, branching is predominantly symmetric (ql ' 0), with little differentia-
tion. This translates to near-exponential growth at the embryonic stage.
• during the bulk of the growth process (6 < l < 30), branching is very close
to balance, although it shifts rapidly from a state of slight imbalance towards
branching (q ≈ 0.45 for l = 6) to a state of slight imbalance towards differentia-
tion (q≈ 0.55 between l = 12 and l = 30). In order to smooth out fluctuations, in
the following, we fit this evolution by a piece-wise function, consisting of a linear
increase in q, up to a constant value (dashed line in Extended Data Fig. 1h).
These results underline a crucial difference between the mammary gland, and
previously studied organs such as the lung, which results in the organ-level hetero-
geneity that one can observe by eye in Fig. 1c [36, 37]. We also found that the
composition of TEBs, defined as the average fraction of EdU+ cells in a EdU+ TEB,
was steady over the entire course of mammary development, lending support to a
model of growth with largely time-invariant properties (Extended Data Fig. 1f).
1.3.2 TEBs of a given level make stochastic decisions as an equipotent pool
To discriminate between the two sub-possibilities described above, we then went into
finer detail and examined whether two offsprings (levels li), sharing a closest ancestor
of level l j, displayed fate correlations. If l j = li− 1, we term the two TEBs “siblings”,
whereas if l j = li− 2, we term the two TEBs “first-cousins”, in analogy with classical
genealogy.
We denote Oi the number of offsprings of a given TEB i (Oi = 0,2). We then
calculated, at each level l, the average number of offsprings 〈O(2)〉(l) of all siblings
(respectively cousins) of TEBs who had two offsprings, as well as the average num-
ber of offsprings 〈O(0)〉(l) of all siblings (respectively cousins) of TEBs who had no
offsprings. If the decision to branch (have two offsprings) or differentiate (have zero
offsprings) is perfectly uncorrelated with previous events (equipotent TEBs), then the
two averages should be the same.
We then calculated the coefficient of fate correlations Fc among either siblings or
cousins at a given level l, defined as
Fc(l) =
〈O(2)〉(l)−〈O(0)〉(l)
2
(5)
where perfect fate correlation translates to Fc(l)= 1, whereas perfect fate anti-correlation
translates to Fc=−1. Indeed, perfect fate correlation entails that the sibling of all TEBs
with 2 offsprings will have 2 offsprings (〈O(2)〉(l) = 2), whereas the sibling of all TEB
with 0 offsprings will have 0 offsprings (〈O(0)〉(l) = 0), leading to Fc(l) = 1.
Notably, we found (Extended Data Fig. 1k) little fate correlation among siblings
(〈Fc(l)〉 = 0.09± 0.19), whereas cousins were slightly anti-correlated, but again with
very small average correlations (〈Fc(l)〉 = −0.09± 0.12). We indicated here respec-
tively the average correlation and standard deviation. This suggests that there are no
preferential regions of growth where TEBs would be predestined to grow more by an
intrinsic bias toward branching.
1.3.3 Predictions on subtree heterogeneity and growth kinetics
If TEB branching follows from such a simple statistical rule, we should be able to
predict the properties of mammary gland topology in full from the average fate proba-
bilities q¯(l) of single TEBs extracted above. Because we identified that the branches
l 6 5 pre-exist at the beginning of puberty, we define a mammary subtree as a branch
s6i of level l = 6 and all of its progeny s
i
j. We then measure the persistence of subtrees
P(l), defined as the probability of finding a subtree reaching a level l, as well as m(l)
the surviving average number of branches in a subtree at a level l.
For a critical birth-death process (q= 1/2), we know that P(l) = 11+0.5l while m(l) =
1+0.5l [41]. Here, we need to account for the fact that the value of q varies according
to the average curve computed earlier. This has been extensively studied as Galton-
Watson processes in varying environments [42]. Therefore, we perform numerical
integration of the modelled dynamics, using the average q(t) measured as an input.
Qualitatively, the initial imbalance towards branching (q< 1/2), and subsequent imbal-
ance towards differentiation causes a correction to the classical Galton-Watson result
(Fig. 1e-f), by creating more mid-sized trees than that expected for neutral dynamics.
Details on the statistical approach used To challenge quantitatively the goodness
of the fit, as well as its explanatory power, we verified the model predictions by calcu-
lating in both cases the coefficient of determination R2, which is the simplest indicator
of the goodness of a fit. In addition, we calculate the standard error of the fit S, which
is an absolute measure of the residuals of the fit, and has been shown to be more
adapted for non-linear fitting procedures [43]. For a good fit, R2 should be as close
to 1 as possible, whereas S should be as close to 0 as possible. Specifically, defining
n as the number of the points being fitted, y¯ = 1n∑
n
i y
i
obs the average of the observ-
able, Stot = ∑ni (yiobs− y¯)2, and Sres = ∑ni (yiobs− yimodel)2, the coefficient of determination
is defined as R2 = 1− SresStot , while the standard error of the fit is given by S=
√
Sres
n .
Model validation Comparing the theoretical persistence with the data, we find ex-
cellent agreement (Fig. 1f), as measured by a coefficient of determination of R2= 0.98,
or alternatively, a standard error of the fit of S = 0.025. Alternatively, given the expo-
nentially decaying nature of the data, we took the logarithms of both the data and the
fit, and again computed a coefficient of determination which puts much more empha-
sis on the tail of the distribution, and still found R2 = 0.89.
As a more intuitive measure of the heterogeneity of subtrees, we then measure
the distribution of total size of subtrees Ti, defined as the total number of branches
having s6i as an ancestor in a given subtree. We then compare this distribution to
the theoretical one (Fig. 1e), and find again a very good agreement between data
and theory, in the absence of any adjustable parameter (as measured by a coefficient
of determination of R = 0.98, or alternatively, a standard error of the fit of S = 0.013).
Again, given the decaying nature of the data, we took the logarithms of both the data
and the fit, and again computed a coefficient of determination which puts much more
emphasis on the tail of the distribution, and found R2 = 0.54.
This suggests that almost all of the size heterogeneity of the mammary gland can
be explained by an extremely simple model of equipotent TEBs, making stochastic
choices in a nearly balanced manner. We therefore do not expect the intrinsic long-
term TEB heterogeneity to be large, although some short-term heterogeneity and
priming could still exist to explain the tail of the size distribution, which is slightly
larger than expected.
Moreover, we can calculate from the stochastic simulation how the fraction of
EdU+ TEBs fEdU changes as a function of TEB level l. Indeed, the fraction of TEBs
that will branch again at a given level is S(l)(1− q(l)), whereas the total number of
existing TEBs is ∑l0 S(l′)q(l′), thus:
fEdU =
S(l)(1−q(l))
∑l0 S(l′)q(l′)
(6)
This assumes that all branches of a given level are generated at a similar time. In
practice, to link this to biological data, we must calculate the relationship between
“real” time and level number. As pubertal mammary development occurs between
3 weeks and 8 weeks, and that up to 25-30 branch levels are generated, we can
approximate that the time between two branching events is approximately 1 days. In
fact, we turned to a 5-week intermediary data point to check this value, and noticed
that gland branch numbers peaked around 15th level at that time, and reaching up to
20th level , thus a growth of around 10 branches in 2 weeks, validating our assumption.
We therefore converted level number into time, and saw that our model predicted
the EdU kinetics (Extended Data Fig. 1g) extremely well (R = 0.96, S = 0.15). This
again strengthens our approach, and shows its predictive value.
In summary, these findings show that the morphology of the mammary gland re-
sults from stochastic fate choices made at the mesoscopic level by TEBs. We demon-
strate indeed that the statistical behaviour of a single TEB can define the organisation
of the full organ.
1.4 Timing of the branching process
In order to go beyond the topology of the mammary gland, one can ask about how the
timing of a branching event is specified. As we cannot access directly the underlying
distribution of timings, we use as a proxy the distance between two branch points,
which reflects the number of cells that were generated between two branch points.
Plotting such distributions are useful, as different modes of regulations of branch-
ing [44] are expected to produce different classes of distributions. For instance, it was
proposed that the lung branching morphogenesis arose from a Turing-like instability
[45], based on the known key signalling factors SHH and FGF10 [46]. This gives rise
to a peaked, Gaussian-like distributions of branch lengths, centered around the Turing
wavelength.
On the other hand, an even simpler model is to assume branching can occur
randomly at any time with equal probability λ , independent of the branch length or on
the history of the system. This is analogous to the well-studied case of radioactive
decay, and yields an exponential branch length distribution:
P(L) = λe−λL (7)
characterised by only one parameter λ , which can be seen as the rate limiting step of
the branching process.
Before asking for the detailed distribution profile, a first important control is to
check whether the average geometrical properties of a duct change throughout de-
velopment. As we note in the main text, we see that both the average duct width
(Extended Data Fig. 1i) and length (Extended Data Fig. 1j) are surprisingly inde-
pendent on the duct level l, as soon as l > 6, i.e. for the part of the tree that is of
interest and was generated during puberty. This echoes the results from the previous
section on the tree topology to argue that pubertal mammary development is largely
a time-invariant process, where the behaviour of a TEB remains largely unchanged
until the end of expansion.
We then examine distributions of branch lengths (Extended Data Fig. 1e), aver-
aging all glands (n = 10 glands, N = 6040 branches). Strikingly, the distribution was
very well described by a single exponential (coefficient of determination R2 = 0.98 and
standard error to the fit S = 0.002), fitting the model of stochastic branching events.
Alternatively, given the exponentially decaying nature of the data, we took the loga-
rithms of both the data and the fit, and again computed a coefficient of determination,
and still found R2 = 0.98, indicative of excellent fits to an exponential. Moreover, the
slope 1/λ of the exponential was highly conserved among different glands, meaning
that all distribution probabilities could be collapsed on a single master curve. This
is particularly surprising as different glands had widely different macroscopic sizes
(2-fold different in total number of branches). However, these differences were not re-
flected at the individual branch level, as we could see only weak correlation between
the total size Tg of a gland g and the slope of its area distribution λg (R2 = 0.3 from an
affine fit).
1.5 Relative contributions of tip branching vs. side branching
Before continuing our analysis, we address the assumption we have made in the main
text, and so far in the Supplementary Note that the bulk of the branching events arise
from a TEB bifurcating into two equal offsprings. Indeed, another source of branching
could co-exist in this system, namely side-branching, that has been described exten-
sively in the context of lung morphogenesis for instance [38]. During side-branching,
an existing duct protrudes outwards to give rise to a new TEB. In mammary morpho-
genesis, the relative amount of side-branching and tip branching has been the subject
of intense scrutiny and debate [47, 48, 49], although the degree of side-branching may
well be dependent on many factors including mouse strain [37]. It is considered that
pubertal branching morphogenesis occurs mainly via tip branching, with a small con-
tribution of side-branching, but that side-branching (so called tertiary as opposed to
the primary and secondary structures deposited during embryogenesis and puberty)
plays a crucial role in young adults, and is under the control of the progesterone-Wnt
pathway [49, 50], responding in particular to the ovarian cycle. Our pubertal findings
are largely in agreement with this picture. As described above and in the main text,
EdU assays show that the large majority of dividing cells reside in TEBs, which dis-
play rapid bifurcations (Extended Data Fig. 1a-c). Ducts, on the other hand, appear
largely quiescent and stationary from our data. This was confirmed by the intra-vital
live imaging experiments, which confirmed the picture from the Edu assays of static
non-proliferative ducts vs dynamic and proliferative TEBs (Fig. 2d).
To challenge this picture further and more quantitatively, we made use of our in-
termediary time point, where we examined mice at 5w of age. We reasoned that if
extensive side-branching occurs between 5w and 8w, in the regions already formed
at 5w (i.e. levels less than 15), the average length between two bifurcation events
should diminish over time. However, comparing the 5w and 8w data, we could find no
statistically significant average branch size differences. Moreover, as we have seen
that branch sizes at 8w obeys a strikingly exponential distributions, we sought to test
whether this was again true at 5w. To assess statistical significance, in this case and
in the rest of the Supplementary Note, we use the Kolmogorov-Smirnov test, which is
a non-parametric test used to compare two sample distributions. It relies on looking at
the maximal distance between two cumulative clone size distribution in two regions.
With f8(x) representing the cumulative distribution of branch areas x at 8 weeks, com-
prising N branches, and f5(x) representing the cumulative distribution of branch areas
x at 5 weeks, comprising N′ branches, we define the Kolmogorov-Smirnov distance
DN,N′(t) = supx| f8(x)− f5(x)|. The two distributions are considered different with a
level of confidence α if (
NN′
N+N′
)
DN,N′(t)> c(α) (8)
where c(α) is the inverse Kolmogorov distribution, from which a P-value can be com-
puted. In comparing the branch length distributions, we had N = 6040 and N′ = 1116,
therefore providing large power to the statistical test. However, we found P > 0.2,
indicative of no statistical differences between 5w and 8w, even at 20% significance
level. This strongly argues against a significant contribution of side-branching to pu-
bertal growth in our mouse strain, although some events could still occur at the very
latest stages of puberty, for the distal-most part of the gland.
1.6 Comparing branching morphogenesis in 4th vs 5th mammary
glands
So far in the analysis, we have grouped together 4th and 5th mammary glands in the
analysis. In principle, their principle designs and architecture could be entirely differ-
ent, and we need therefore to check a posteriori that such a grouping made sense
from a statistical perspective. As shown on Extended Data Fig. 3d-g, we grouped sep-
arately the 4th and 5th glands (n=5 in both cases), and noticed that both types display
time-invariant average branch length and width (Extended Data Fig. 3d). Moreover,
examining the average probability of TEB termination at a given segment level (Ex-
tended Data Fig. 3e), revealed that both gland types display similar convergence to a
near-balance between termination and branching. In order to probe the similarity be-
tween the two types more quantitatively, we then computed the subtree persistence
and size distribution for each, and calculated the Kolmogorov-Smirnov distance for
each quantity (n= 86 subtrees for 4th glands and n= 74 subtrees for 5th glands). For
both the subtree persistence (Extended Data Fig. 3f) and size distributions (Extended
Data Fig. 3g), we found no statistically significant differences (P> 0.2 in both cases)
even at 20% significance level, validating our approach of grouping both glands to-
gether in the analysis.
1.7 Comparing branching morphogenesis at various Tamoxifen
level
Moreover, we have also grouped together mammary glands that were untraced (no
tamoxifen injections) and mammary glands which were injected with various doses of
Tamoxifen (0.2 and 1.5 mg). As detailed in the main text and shown on Extended Data
Fig.3a, we have optimised the dose of tamoxifen so that we can get confetti clones
without interfering with the branching pattern. We have determined that a dose of 3mg
of Tamoxifen caused a delay in branching morphogenesis, so we restricted ourselves
to lower doses. In particular, we show in Extended Data Fig. 3a-c that low doses
of Tamoxifen (0.2 and 1.5 mg) did not impact on the time-invariance of the branch
length and width (Extended Data Fig. 3b), nor did it change the growth characteristic
q(l) of the mammary gland (Extended Data Fig. 3c). Once ungrouped for Tamoxifen
doses, the number of subtrees in each case was relatively small, diminishing slightly
the statistical power of the test (n = 17 for 2 mice for un-injected, n = 26 for 2 mice
for 0.2mg and n = 68 for 5 mice for 1.5mg). Overall, we still compared 0 and 0.2mg
to 1.5mg, both relative to their subtree size distribution and subtree persistence. In
the 0 vs 1.5mg comparison, we found respectively P > 0.2 and P > 0.2. In the 0.2 vs
1.5mg comparison, we found respectively P > 0.1 and P > 0.05. Overall, this shows
that Tamoxifen does not have a discernible effect on the branching morphogenesis at
these doses.
1.8 Applicability of these findings to kidney morphogenesis
We then sought to determine whether some similar design principles could in principle
explain the ductal morphogenesis of other organs such as kidney. Although kidney
was thought to be a rather stereotypical organ [51], recent findings have suggested
an un-appreciated level of complexity by performing detailed three-dimensional recon-
structions [52]. As detailed in the Materials and Methods, we thus performed a similar
branching analysis on an E16 murine kidney (Extended Data Fig. 10a,b). Interest-
ingly, we first noted that although kidney undergoes a longer period of pure symmetric
branching, the probability for branches to terminate q(l) still displays a convergence
towards a state nearly balanced between termination and branching (Extended Data
Fig. 10c), very reminiscent of the paradigm we uncovered in mammary gland, and in
agreement with the findings of Ref. [52].
Before proceeding to simulate tree heterogeneity, we first sought to check whether
the branch length distribution in kidney was an exponential, indicative of a stochastic
branching process. Interestingly, we could indeed observe that the distribution was
extremely well fitted by a single exponential (Extended Data Fig. 10d), again mirroring
our findings in mammary gland (one should note here that we excluded from the anal-
ysis the first five levels, which were on average longer than the other segments). We
then calculated the coefficient of fate correlations Fc among either siblings at a given
level l, using the same definition as before, i.e. where perfect fate correlation trans-
lates to Fc(l) = 1, and perfect fate anticorrelation translates to Fc = −1. Interestingly,
we found 〈Fc(l)〉= 0.28±0.27) (average correlation and standard deviation), indicative
of existent, albeit relatively weak and noisy fate correlations. One should note that the
existence of such a non-zero correlation could arise from the intrinsic anisotropy of
growth present in the kidney, where different axes branch and elongation at different
rates [52].
We next performed a Galton-Watson modelling of subtree heterogeneity, assum-
ing as in mammary that tips were all equipotent, and branched/terminated stochasti-
cally, according to the average properties measured in Extended Data Fig. 10c (we
defined a subtree as all offsprings of branches of level 7). Strikingly, we found that
the subtree heterogeneity, both measured from the subtree persistence (Extended
Data Fig. 10e) and size distribution (Extended Data Fig. 10f), fitted extremely well
to the observed data. Both experimental quantities displays broad exponential tails,
again reminiscent of mammary structural heterogeneity. This confirms that although
there might be some bias linked to the directionality of growth, a model of equipotent
tips making stochastic fate choices still account well for the subtree heterogeneity
displayed in kidney.
1.9 Discussion and outlook on the balance between branching
and termination
Taken altogether, our findings above provide a quantitative basis for the often-made
qualitative statement that mammary gland structure is disorganised and not stereo-
typical [36, 53] compared to other branched organs. A key question that emerges
from this observation is therefore the mode of regulation of the branching probability
q(l), and whether the value of 1/2, characteristic of organ level balance represent a
bona fide attractor of the system. Although past literature has identified a multitude
of factors that affect branching [36, 37, 53, 54] (including hormones, stromal factors
such as MMPs, and epithelial signalling pathways such as FGF), the fact that q(l) is
positioned so close to balance for such a long time, as well as the time-invariance of
the branching process and geometry in general, suggests the existence of a simple
underlying rule governing this dynamics.
One way to explain such dynamics from a theoretical perspective would be to con-
sider branching as a (stochastic) default state of a TEB, with associated rate constant
b0, and termination as resulting from a crowding effect, dependent on the density
of TEBs at a given time t, or segment level l, which we note n(l). As discussed
above, this crowding effect could theoretically arise from direct TEB-TEB interactions,
or from competition for resources. One can then write the mean field equation for
the evolution of the density n of TEBs as n′(l) = n(l)(b0− f (n(l)) where f is a generic
monotonously increasing function which characterises the dependency of termination
upon TEB density. This is a logistic equation which, assuming that liml→∞ f (l) > b0,
predicts a stable steady state for n = n0. Interestingly, a constant density n(l) trans-
lates into q= 1/2, showing that such generic “crowding“ regulation of TEB termination
could explain that the balance between branching and termination is an attractor of
branching morphogenesis.
However, the equipotency of TEBs on at organ level that we have uncovered in
this section still leaves open the key question of the potency and hierarchy of individ-
ual MaSCs in a given TEBs. In principle indeed, each TEB could be organised either
in a hierarchy with a very small number of true MaSCs giving rise to a large pool of
committed progenitors, either as an equipotent population of cells with equal popula-
tion for morphogenesis. We therefore now turn to the clonal dynamics at the cellular
levels which are produced by the mesoscopic branching choices made by TEBs.
2 Clonal dynamics
As stated in the main text, the time-invariance of the properties of the TEBs and ducts
argue for a simple model of branching at the cellular scale. Based on the clonal data,
in the following, we consider independently the basal and luminal compartments. We
therefore assume that tips contain N (basal or luminal) stem cells, which segregate
into two approximately equal parts (N/2) at each budding event, before expanding
through symmetrical division to reconstitute the N stem cells that characterise the
TEB after budding.
Upon confetti-lineage tracing, n stem cells are labelled in a given colour in the TEB.
But how are these cells repartitioned in the two daughter TEBs following bifurcation,
assuming that where a stem cell will go depends only on its location at the time of
budding, and not on its history (a Markovian process). More generally, since stem
cells can migrate within the TEB between two consecutive branching events, leading
to the fragmentation and merging of clones, a generic rule on clonal evolution cannot
be derived. However, if we define D as the effective diffusion coefficient of a stem
cell in the TEB, L≈ 200µm the TEB size along the axis of the duct and τ ≈ 2 days the
characteristic time between consecutive branching events, analytical expressions can
be obtained for the two limiting cases of low diffusion (no-mixing) D L2/τ and high
diffusion D L2/τ (well-mixed). We therefore successively examine the two cases:
2.1 Well-mixed model
In this regime, we assume that the diffusion constant is sufficiently large that the po-
sition of stem cells between two consecutive budding events are uncorrelated. This
means effectively that the N/2 stem cells picked to go in one of the daughter TEBs
are randomly chosen within the pool of N stem cells. This is a version of the classic
“occupancy problem”, in the case of a drawing without replacement and with the con-
straint that each tip has to have N/2 stem cells. In this case, one can show that the
probability Πn→k to obtain k stem cells in a daughter TEB, given that one had n stem
cells in the mother TEB is:
Πn→k =
(
n
k
)
(N−n)!
N!
(N/2)!
(N/2− k)!
(N/2)!
(N/2−n+ k)! (9)
In the limit of an infinitely large stem cell compartment (N→ ∞), this formula reduces
to the familiar binomial distribution:
Πn→k =
1
2n
(
n
k
)
(10)
One can then write the Master equation associated with the above dynamics. We
define the probability Pn(t) of having, at the generation t, n stem cells in a TEB. At the
next budding (generation t+1):
Pn(t+1) =
N
∑
i=n/2
Πi→n/2Pi(t) (11)
The sum starts at i= n/2 because lower stem cell numbers than n/2 cannot contribute
to Pn(t + 1). Moreover, the probability Πi→n/2 considers selecting n/2 cells out of i,
because of the subsequent doubling in numbers, which we assume instantaneous
here, and will therefore result in n labelled cells after bifurcation. The initial condition
is Pn(0) = δn,ni, i.e. ni stem cells of a given colour is clonally induced at the beginning
of the process. We first start with the simplest case ni = 1. As the overall mean clone
size is constant (〈P〉(t) = ∑∞i=0 iPi(t) = ni), a more interesting quantity is the surviving
mean clone size m(t) defined as
m(t) =
∑∞i=0 iPi(t)
1−P0(t) =
1
1−P0(t) (12)
2.1.1 Case of very large number of stem cells N→ ∞
To solve the Master equation, one may note that the time evolution for P0(t) is simply
given by
P0(t+1) = P0(t)+
∞
∑
i=1
Pi(t)
2i
(13)
The dynamics is analogous again to a Galton-Watson process [55] with, from the
perspective of a given offspring TEB, each cell choosing independently and stochas-
tically to either duplicate (i.e. go in that TEB) or become extinct (go into the other
offspring TEB), with probabilities 1/2 and 1/2. With τ the time between two budding
events, the mean persisting clone size, then converges towards m(t)→ t2τ . Extended
Data Fig. 8b shows a representative example of such dynamics for three initial con-
ditions (ni = 1,2,3). One can indeed verify that the curves converge rapidly to the
predicted linear slope.
2.1.2 General case of finite N
For a finite number of stem cells, the above expression is still valid at short times,
defined as m(t) N. At longer times, one gets monoclonal conversion of TEBs,
and one can get analytical limits for N −m(t) N by applying precisely the same
expression as above for the non-labelled cell fraction of persistent size N−m(t).
We now turn to the no-mixing model to contrast the results we have obtained.
2.2 Low-mixing, spatially structured model
In the low mixing regime, no clonal fragmentation occurs, so that spatial correlations
have to be taken into account. We make the same assumption of n labelled stem
cells in a total N, but we now take into account spatial structure, so that each stem
cell occupies a position i∈ [1,N] on a torus. We define an indicator function ni(t), such
that ni = 1 if a cell is labelled, and ni = 0 otherwise. When a budding event occurs,
N/2 contiguous cells are chosen, starting at a random position j, and duplicated, so
that {
n2k+1(t+1) = n j+k(t)
n2k+2(t+1) = n j+k(t)
for k ∈ [0,N/2−1]. (14)
We start again with the initial condition ni(0) = δi,1. As stated above, a fundamental
property of the system is that labelled clones then remain cohesive, and that the mean
number of cells in a TEB m = ∑Ni=1 ni is even. Working out again the combinatorics,
one can arrive to simple expression for the transition rates Πn→k under the condition
k < N/2, depending on the value of k relative to n:
Πn→2k>2n = 0
Πn→2n =Πn→0 =
N/2−n+1
N
Πn→2k<2n = 2N
(15)
The first line stems from the obvious fact that stem cell number cannot increase by
more than a factor 2 per budding. The second line deals with the case when all
labelled cells are placed in the same TEB. The third line deals with boundary events,
which can only occur 2 times out of N, i.e. when the selected position is k positions
away from the boundary of the clone. The Master equation for Pn(t) then be written{
P0(t+1) = P0(t)+∑
N/2
k=1
N/2−k+1
N P2k(t)
Pm(t+1) =
N/2−m/2+1
N Pm/2(t)+
2
N ∑
N/2
k=m/2+1P2k(t)
(16)
so that the evolution of P0(t) is dictated by a simple equation, as in the continuous
limit:
τ
dP0
dt
=
N/2+1
N
(1−P0(t))− 12N (17)
2.2.1 Case of very large number of stem cells N 1
For N → ∞, one expects an exponential clonal growth, as the only possible choices
are either not selecting any labelled cells (probability 1/2), or selecting all of them
(probability 1/2). The mean surviving clonal size and clonal persistence in the scaling
regime are indeed
m(t) = ekt (18)
and
P0(t) = 1− e−kt (19)
where k = ln(2)/τ is the rate of TEB bifurcation. This is fundamentally different from
the well-mixed model, which displayed a linear growth regime for N→∞. As expected,
the persistence converges at long times to 1−P0(t)→ 1/N, as the experiment marked
initially 1 cell out of N.
2.2.2 General case of finite N
For finite N, we need to also calculate the transition probabilities for n > N/2. The
same combinatorics argument as before can be used to calculate the transition prob-
abilities: 
Πn→2k<2n−N = 0
Πn→2k=2n−N =
n−N/2+1
N
Πn→2k>2n−N = 2N
Πn→N = n−N/2+1N
(20)
We thus approach monoclonal conversion (N −m(t)→ 0) exponentially. Extended
Data Fig. 8a,b (center) shows a representative example of tree reconstruction with
such dynamics, for N = 200. We also plot the evolution of the mean persisting clone
size in time (Extended Data Fig. 8c) to verify that we indeed get two successive
phases of exponential growth and plateauing.
2.3 Data analysis
Although the theoretical modelling above provides predictions for the evolution in time
(or segment level) of clonal composition within TEBs, the corresponding data is not
directly accessible experimentally, as one would need to sacrifice mice nearly every
day to monitor the evolution of the TEB composition. However, as discussed in the
main text, the ducts provide a historical record of the face choices made by TEBs at
time time that the duct was formed. Therefore, by combining tree reconstruction and
spatial positioning of all labelled cells within the tree, we can infer the composition
of TEBs at all past time points. Two crucial hypotheses underlying this assumptions
are that the morphogenesis is tip-driven and that ducts are static once formed. As
discussed in the main text and in the sections above, we verified the validity of these
assumptions via multiple independent experiments: Edu assays (Fig. 1a and Ex-
tended Data Fig. 1), continuous intravital live imaging of TEBs and ducts (Fig. 2d) as
well as comparison between the branch structure and geometry in 5 vs 8 week-old
mice (Extended Data Fig. 1e).
2.3.1 Calculation of MaSC number
As the MaSC number N in a given TEB is an important input of the modelling, we
briefly summarise how we obtained it experimentally. We only made use here of
the ultra-low dose (0.2mg Tamoxifen) lineage tracing experiments from 3 to 8 weeks
of age, where clonality is insured. We then defined as before a subtree as all the
branches originating from a branch of level 6. In all subtrees considered, we either
found not a single confetti labelled cell, or multiple confetti labelled cells throughout
the subtree, all sharing the same colour and cell type (basal or luminal). In total, we
found 19 subtrees which were labelled (11 with luminal cells and 8 with basal cells).
For each cell type, we then calculated the average labelled cell fraction in a subtree,
i.e. the number of labelled cells divided by the total number of cells in that subtree. We
counted manually labelled cells, and used the sum of all branch area as a proxy for
total cell number, having measured independently the average cellular density of both
basal and luminal cells (Extended Data Fig. 2h). Since we know that all of the labelled
cell in a given subtree originate from a single cell, and we are tracing representatively
the whole population, we can deduce the effect number of MaSCs in a TEB by taking
the inverse of the labelled cell fraction (Extended Data Fig. 6a for a schematic). In
the main text, we thus report the mean and s.e.m. of the number of basal and luminal
MaSCs (Extended Data Fig. 6d). Notably, both numbers were very large, indicating
that the limit of N→ ∞ in the previous section is a relevant one. Moreover, this large
number provides an explanation as to why we never see monoclonal ducts, as the
system does not have enough time to reach monoclonality. This leads us, in the
following sections, to consider in more detail the quantitative cellular dynamics during
branching and elongation.
2.3.2 Tracing from onset of puberty to 8 weeks
Firstly, in contrast to the theory, our full-gland reconstruction does not give the clonal
composition of all TEBs at every time point of development, but rather the ductal
imprint left by these TEB clones. Moreover, because the length and width of a given
branch is stochastic, and can vary by a large amount, the relevant proxy to use is the
clonal density ρi(l) of a given colour in a branch si of level l, rather than the absolute
number of cells labelled in that branch.
We start by considering the 3 to 8 week tracing (Fig. 4 and Extended Data Fig.
2a-c). We calculate for all glands of low- and medium-density labelling the average
clonal density ρ(l) of either basal and luminal cells at a level l, and normalise it by
their respective values for l = 6, since we assume all initial labelling events occurred
at this stage. Summing all clones, we had a total of 157 clones of basal cells and
449 clones of luminal cells to analyse, ranging from level numbers between l = 6 and
l = 30.
From the resulting graph (plotted on Fig. 4a), it is apparent that the increase in
mean clone size as a function of branch number is very slow, fitting with a linear
rather than exponential increase. Moreover, although the dynamics are rather noisy,
the slope of the increase is close to the value predicted in the well-mixed limit, show-
ing that the mammary gland development operates in that regime. This provides a
quantitative basis for our observation of a neutral drift from the proximal to distal end
of the tree (Fig. 4 and Extended Data Fig. 4,9)
As we have seen before, the one fitting parameter left in order to predict quantita-
tively the full curve is the initial condition, i.e. the number of labelled cells at the first
branching event. Indeed, we have seen that the average timing between branching
events is around T ≈ 1.5 days, so that on average, we clonally label a cell 1 day before
the first branching event of its TEB. As we have inferred a division rate around 12h,
we would therefore expect to have 2-4 labelled cells initially. Alternatively, a possi-
ble source of variation in this number is the initial condition of the dynamics. In our
idealised model, we have assumed that all active TEBs are all initially at the same
level. In fact, as supported by the Edu assays, we see that the segment level of ac-
tive TEBs is distributed around 6, but with a significant amount of variation. This will
result in some TEBs reaching level 6 after having already performed several rounds
of bifurcation, and expanded their clone sizes.
Interestingly, when determining the optimal fit of the data, we found an initial con-
dition of n= 3 cells (yielding a coefficient of determination of R2 = 0.65 and a standard
error to the fit S = 0.69). Moreover, as an important consistency check, we have also
calculated the labelled cell fraction as a function of segment level, which is the product
of the mean clone size calculated above, and the clonal persistence (i.e. the fraction
of branches of a given level which contain at least one confetti labelled cell). If the
cells targeted by the lineage tracing really function as an equipotent population which
is representative of the tissue, then one expects the labelled cell fraction to remain
to its initial value for later segment levels. Formulated alternatively, given that mean
clone size increases linearly, this predicts that persistence should decrease in inverse
proportion, as in the classical critical Galton-Watson case. We therefore turned to the
data and show (Extended Data Fig. 2c) that the labelled cell fraction indeed hovers
around its initial value as segment level 6, consistent with the representative labelling
of an equipotent population undergoing neutral drift due to bifurcation events.
Moreover, the theory predicts that the clone size distributions Pn(l) should con-
verge rapidly to single-exponentials, reflecting the equipotency of MaSCs in a TEB.
This scaling behaviour is expressed as
Pn(l) =
1
〈n(l)〉 f
(
n
〈n(l)〉
)
where 〈n(l)〉 is the average clone size at a level l, and f (x) = exp(−x). One should
note that the convergence to this scaling form is not affected by the aforementioned
fact that the initial condition may vary in different TEBs. Because each level contains a
relatively small number of clones, it was impossible to test whether each level showed
exponential distributions. Therefore, we use a 3-level binning to achieve better statis-
tical significance, and display on Fig. 3b-c the distributions, rescaled by their average
clone size. We find consistently that the distributions display scaling behaviour and
are well-fitted by a single exponential, as expected from an equipotent population. To
make this analysis more quantitative, we again assess P-values for the exponentiality
of the distributions at all segment numbers, using as before a Kolmogorov-Smirnov
test. One should note that the classical Kolmogorov-Smirnov test assumes statistical
independence between the model and the experimental data. Here, we have used
some information on the mean clone size distribution, and simply want to know if the
rescaled form of the distribution conforms to an exponential. Therefore, we made use
here the modified tables of critical values of the statistic Dn, testing whether the ex-
perimental data follows an exponential distribution when the mean is estimated from
the data (following the classical results of Lilliefors [56] and others).
For the basal distributions, at all 7 segment levels plotted in Fig. 3b, the P-values
are respectively:
P= [> 0.2,0.05,> 0.2,> 0.2,> 0.2,0.04,> 0.2]
For the luminal distributions, at all 7 segment levels plotted in Fig. 3b, the P-values
are respectively:
P= [> 0.2,> 0.2,> 0.2,> 0.2,> 0.2,> 0.2,> 0.2]
This shows that the rescaled distributions are consistently well-fitted by a single
exponential at all time points. The combination of exponential clone size distributions,
and linear mean clone size increase argue strongly in favour of MaSCs of each lin-
eage working as a single well-mixed compartment. One should note however that due
to the large size of the compartment, neutral drift is not sufficient to ever reach TEB
monoclonality in the time frame of pubertal expansion, making it therefore impossi-
ble to assess the size of these compartments from the neutral drift and distributions
alone (as is done typically for neutral drift in intestinal stem cells). It could therefore
be possible that some large independent compartments exist within a TEB lineage
(i.e. within the constraint of basal and luminal unipotency), although they would each
need to be large, have similar sizes, as well as similar growth and mixing dynamics in
order to fit our data.
2.3.3 Pubertal tracings at different time points
In order to verify the assumption of our model, and the fact that branching morpho-
genesis occurs in a time-invariant manner, we have performed two other tracings: the
first from 3 to 5 weeks, in order to check that the first half of the tree deposited during
the first half of morphogenesis shows the same clonal property than the complete
trancing, and the second one from 5 to 8 weeks, in order to verify that we still label an
equipotent population in a representative manner at 5 weeks. The first consistency
check we performed is to compute the labelled cell fraction in both cases. Interest-
ingly, and similarly to the complete tracing, we noticed that the rescaled labelled cell
fraction remains on average unchanged relative to its original value, hovering around
1 in all cases (Extended Data Fig. 2d,e). This shows in particular that even at 5
weeks, we still label representatively the MaSC population.
Moreover, we then looked at the 3 to 5w tracing in more details, and found that
both luminal and basal mean clone sizes increases together as a function of segment
level, in a similar amount to the mean clone size increase of the complete 3 to 8w
tracing (Extended Data Fig. 2f from n=3 mammary glands). This again validates
independently our results and modelling approach. In the 5 to 8 w tracing, it was not
possible to plot consistently the normalised mean clone size, because of a technical
issue: at 5w, in contrast to 3w, the initial condition is extremely ill-defined, with active
TEBs ranging from level 8 to level 20. This makes a grouping of subclones size as
a function of segment level ill-defined from a statistical perspective, although the fact
that the labelled cell fraction is constant informs us on the representativeness of the
equipotent population.
We then examined again the rescaled clone size distributions in the 3 to 5 week
tracing (Extended Data Fig. 2g) , both for the basal and luminal cells, to verify whether
these converged to an exponential. We again used to 3-level binning to increase
statistics, and examined the levels 9-11, 12-14, 15-17 and 18-20. As before, we
calculate P-values from the modified Kolmogorov-Smirnov test, in order to examine
whether the distributions can be distinguished from exponentials. For the basal distri-
butions, at the 4 segment levels plotted on Extended Data Fig. 2g, the P-values are
respectively:
P= [> 0.2,0.05,> 0.2,< 0.01]
For the luminal distributions, at the 4 segment levels plotted on Extended Data Fig.
2g, the P-values are respectively:
P= [> 0.15,< 0.01,> 0.15,> 0.1]
Again, in all but two cases, which could reflect statistical or mouse-to-mouse varia-
tions, this shows that the rescaled distributions are well-fitted by a single exponential.
Combined with the 3 to 8 week tracing, this shows that the data conforms well to our
model of a single equipotent population of MaSCs in each TEB, undergoing extensive
mixing and being segregated stochastically.
2.4 Clonal loss
So far, we have assumed that MaSCs in a TEB at generation i are always transmitted
to either one of the two offspring TEB at generation i+ 1. This is clearly an approx-
imation, as one would expect symmetrically dividing MaSCs in a TEB to compete
for space, resulting in an additional source of neutral drift. Note however that such
phenomena is hard to quantitatively distinguish from segregation-induced neutral drift
on the basis of clone size distributions alone, as both cause exponential clone size
distributions and linear surviving clone size increase in time. However, using de-
tailed reconstructions, we can measure at each level the probability for a given colour
and cell type to be found again in at least a segment of the remaining subsequent
branches of the tree. This gives a measure of clonal loss (we exclude in this calcula-
tion very short terminal TEB branches, as it is hard to resolve to which segment a cell
belongs to in this case). Quantitively, we indeed see in the data a non zero, although
low, probability of clone dropping (Extended Data Fig. 4f), which is largely confined
to low-generation numbers, i.e. zones close to the rudimentary tree. This is what
is expected theoretically, as well-mixed clones reaching larger sizes have a vanish-
ingly small probability to drop out, although an exact calculation of the loss probability
would require a more detailed understanding of the exchange rates between different
TEB locations. A technical issue which complexifies the analysis here is that such
early ”dropping out” is very hard to distinguish from having hit TA cells initially. This is
reinforced by the fact, discussed in the sections above, that we have imperfect knowl-
edge of the initial condition of the problem, i.e. the initial generation of each given
TEB upon clonal induction.
2.5 Short-term dynamics
As discussed in the main text, a striking feature of the short-term clonal data is the
fact that only cells located at the border of the TEB contributed to the elongation of
the connected duct, whereas cells located at the tip of the TEB did not (Fig. 3). This
could be consistent with two classes of behaviour. On one hand, tip cells could form
a distinct population from border cells with the former describing stem cells and the
latter operating as transit-amplifying (TA) cells. On the other hand, all cells in the TEB
could function long-term as a single equipotent population, with local spatial cues
providing a short-term positional fate bias. Our data suggests strongly that the latter
option is the most likely.
First of all, we have seen from full clone reconstruction (Fig. 2a-c and Extended
Data Fig. 6), that TEBs must consist of a large number of functional stem cells, match-
ing rather well the sum of border and tip cells (as counted from 3D reconstructions).
Second, labelled confetti MaSCs seem to be randomly distributed over the TEB with-
out any bias towards a certain position, and a large fraction (≈ 50%) of single-confetti
cells in TEBs were located in the border region. If this border region was a TA com-
partment, one would expect a very large rate of clonal loss, which does not occur.
Rather, the ductal record displays “intermittency” (Extended Data Fig. 4e), a fea-
ture we have not analysed yet within our framework. We define intermittency as a
disconnection in the clonal record in the duct - for example, on Fig. 3d red cells in
branches 5-10, no labelled cells in branches 10-12, and red cells again in branches
12-20. Quantitatively, we define Pi(l) the probability of having an intermittent be-
haviour lasting for l branches. Extended Data Fig. 4e shows that this probability
distribution. Strikingly, the distribution of intermittency is rather well-fitted by a single
exponential.
A first very simple hypothesis for this behaviour would be that the number of stem
cells in a TEB is large compared to the number of cells in a typical duct, so that not
every stem cells in a TEB need to divide in order to build a given branch. However, this
does not seem to be the case, as a typical TEB contains around 250 cells, whereas a
branch contains on average 750 cells. Moreover, intermittency does not impact upon
short branches only, implying that it results from the actual behaviour of stem cells,
rather than non-representativeness of the ductal record.
Interestingly, we note that, in order to make the connection with the previous sec-
tion, a large, but not infinite, degree of mixing would yield exactly this short-term
intermittency: tip cells would need a certain time tD ≈ L2/D to contribute to a duct,
and if tD is close to the time between two branching events τ, one could explain this
phenomenon of short-term intermittency, while maintaining the long-term clonal dy-
namics close to the well-mixed model.
References
[36] Gjorevski, N., and Nelson, C. M. (2011). Integrated morphodynamic signalling of
the mammary gland. Nature Reviews Molecular Cell Biology, 12(9), 581-593.
[37] Sternlicht, M. D., Kouros-Mehr, H., Lu, P., and Werb, Z. (2006). Hormonal and
local control of mammary branching morphogenesis. Differentiation, 74(7), 365-
381.
[38] Metzger, R. J., Klein, O. D., Martin, G. R., and Krasnow, M. A. (2008). The
branching programme of mouse lung development. Nature, 453(7196), 745-750.
[39] Simons, B. D., and Clevers, H. (2011). Strategies for homeostatic stem cell self-
renewal in adult tissues. Cell, 145(6), 851-862.
[40] Krieger, T., and Simons, B. D. (2015). Dynamic stem cell heterogeneity. Devel-
opment, 142(8), 1396-1406.
[41] Harris, T. E. (2002). The theory of branching processes. Courier Corporation.
[42] Jagers, P. (1974). Galton-Watson processes in varying environments. Journal of
Applied Probability, 174-178.
[43] Brown, A.M. (2001). A step-by-step guide to non-linear regression analysis of
experimental data using a Microsoft Excel spreadsheet. Computer methods and
programs in biomedicine 65, no. 3: 191-200.
[44] Iber, D., and Menshykau, D. (2013). The control of branching morphogenesis.
Open biology, 3(9), 130088.
[45] Menshykau, D., Kraemer, C., and Iber, D. (2012). Branch mode selection during
early lung development. PLoS Comp Biology, 8(2), e1002377.
[46] Metzger, R. J., and Krasnow, M. A. (1999). Genetic control of branching morpho-
genesis. Science, 284(5420), 1635-1639.
[47] Visvader, J. E., and Stingl, J. (2014). Mammary stem cells and the differentiation
hierarchy: current status and perspectives. Genes and development, 28(11),
1143-1158.
[48] Sternlicht, M. D. (2006). The cues that regulate ductal branching and morpho-
genesis. Breast Cancer Res, 8(1), 201.
[49] Robinson, G. W., Hennighausen, L., and Johnson, P. F. (2000). Side-branching
in the mammary gland: the progesterone?Wnt connection. Genes and develop-
ment, 14(8), 889-894.
[50] Atwood, C. S., Hovey, R. C., Glover, J. P., Chepko, G., Ginsburg, E., Robison,
W. G., and Vonderhaar, B. K. (2000). Progesterone induces side-branching of
the ductal epithelium in the mammary glands of peripubertal mice. Journal of
Endocrinology, 167(1), 39-52.
[51] Al-Awqati, Q., and Goldberg, M. R. (1998). Architectural patterns in branching
morphogenesis in the kidney. Kidney international, 54(6), 1832-1842.
[52] Sampogna, R.V., Schneider, L. and Al-Awqati, Q., 2015. Developmental Pro-
gramming of Branching Morphogenesis in the Kidney. Journal of the American
Society of Nephrology, pp.ASN-2014090886.
[53] Ochoa-Espinosa, A., and Affolter, M. (2012). Branching morphogenesis: from
cells to organs and back. Cold Spring Harbor perspectives in biology, 4(10),
a008243.
[54] Inman, J. L., Robertson, C., Mott, J. D., and Bissell, M. J. (2015). Mammary
gland development: cell fate specification, stem cells and the microenvironment.
Development, 142(6), 1028-1042.
[55] Athreya, K. B., and Ney, P. E. (2012). Branching processes (Vol. 196). Springer
Science and Business Media.
[56] Lilliefors, H.W. (1969). On the Kolmogorov-Smirnov test for the exponential dis-
tribution with mean unknown. Journal of the American Statistical Association 64,
no. 325: 387-389.