ABSTRACT
Branching morphogenesis is a characteristic feature of many essential organs, such as the lung and kidney, and most glands, and is the net result of two tissue behaviors: branch point initiation and elongation. Each branched organ has a distinct architecture customized to its physiological function, but how patterning occurs in these ramified tubular structures is a fundamental problem of development. Here, we use quantitative 3D morphometrics, time-lapse imaging, manipulation of ex vivo cultured mouse embryonic organs and mice deficient in the planar cell polarity component Vangl2 to address this question in the developing mammary gland. Our results show that the embryonic epithelial trees are highly complex in topology owing to the flexible use of two distinct modes of branch point initiation: lateral branching and tip bifurcation. This non-stereotypy was contrasted by the remarkably constant average branch frequency, indicating a ductal growth invariant, yet stochastic, propensity to branch. The probability of branching was malleable and could be tuned by manipulating the Fgf10 and Tgfβ1 pathways. Finally, our in vivo data and ex vivo time-lapse imaging suggest the involvement of tissue rearrangements in mammary branch elongation.
INTRODUCTION
Branching morphogenesis is a fundamental developmental process driving the formation of a number of organs and tissues, such as the mammary gland, salivary gland, kidney, lung, blood vessels and nervous system (Goodwin and Nelson, 2020; Lang et al., 2021; Lu and Werb, 2008). Branching involves repeated events of branch point formation and tube elongation, which enable the maximization of the functional area of an organ in a limited three-dimensional (3D) space. Two distinct mechanisms of branch point formation have been described: tip splitting (or clefting) and side (lateral) branching (Lu and Werb, 2008; Myllymäki and Mikkola, 2019). Each organ has its unique way of generating the branching pattern – choice of branch site initiation, and the length, diameter and spacing of branches – that gives rise to an organ-typical structure adapted for its specific function. For example, lung branching morphogenesis is stereotyped. It involves three geometrically distinct modes of branching: domain (lateral) branching, planar bifurcation and orthogonal bifurcation. Each of these modes occurs in a fixed order; thus, lung branching is strictly regulated via pre-determined sequential branching events (Metzger et al., 2008). Similar to lungs, kidney branching morphogenesis is highly reproducible (Short and Smyth, 2016). The ureteric bud undergoes a series of reiterative rounds of terminal bifurcations early and trifurcations later, resulting in a ureteric tree where the extent and gross pattern of branch elaboration is conserved (Short et al., 2014). Salivary gland morphogenesis also proceeds without side-branching events and instead involves repetitive clefting of existing buds that widen rather than elongate before the next round of clefting events (Hauser and Hoffman, 2015).
In contrast to many other organs, the branching of the mammary gland appears stochastic (Macias and Hinck, 2012; Scheele et al., 2017). It is also unique in that it generates an ‘open’ architecture, which is a relatively sparse ductal network that leaves space for pregnancy-induced milk-producing alveoli. How the macroscopic features of the mammary gland, including the size, network topology and spatial patterning, are designated remains poorly understood. Branching of the mammary gland takes place in distinct developmental stages, each with discernible features: embryonic/pre-pubertal, pubertal and reproductive (Macias and Hinck, 2012). In mice, five pairs of small ductal trees with 10-20 branches form before birth by unknown mechanisms (Spina and Cowin, 2021). Initially, mammary rudiments branch as a solid mass of cells that, soon after birth, resolve into a bilayer of inner luminal and outer basal cells. Thereafter, the ductal system grows isometrically until the onset of puberty, when ductal tips turn into larger multi-layered bulbous terminal end buds (TEBs) that orchestrate the ensuing burst of ductal growth (Paine and Lewis, 2017). Two-dimensional morphometric studies of fixed samples, together with mathematical modeling, suggest that branching during puberty is mainly driven by TEB bifurcations, whereas lateral branches emerge from previously formed ducts during the estrous cycle and pregnancy (Brisken and Ataca, 2015; Scheele et al., 2017).
Branch pattern formation is governed by reciprocal tissue crosstalk between the epithelial tree and the surrounding mesenchyme, mediated by both shared and organ-specific secreted signaling molecules and extracellular matrix (ECM) components (Goodwin and Nelson, 2020; Lang et al., 2021). Embryonic mammary gland branching occurs in a mesenchymal fat pad precursor tissue that matures postnatally into a complex stroma consisting of adipocytes, fibroblasts, immune and endothelial cells – all with established roles in morphogenesis (Paine and Lewis, 2017; Wang et al., 2020). Mammary gland development during puberty and reproductive life crucially depends on systemic hormonal cues, the effects of which are, however, mediated by locally produced paracrine and autocrine factors, thereby paralleling embryonic branching morphogenesis (Brisken and Ataca, 2015; Macias and Hinck, 2012). Many of these factors have been recognized, including fibroblast growth factors (Fgf), Wnt ligands, transforming growth factor β1 (Tgfβ1), and epidermal growth factor (Egf) and tumor necrosis factor (Tnf) family members (Ciarloni et al., 2007; Luetteke et al., 1999; Macias and Hinck, 2012; Paine and Lewis, 2017; Voutilainen et al., 2012). However, what is less well understood is to what extent these pathways regulate growth and/or differentiation versus branch patterning per se. For example, although epithelial Fgf receptor signaling is pivotal for ductal branching, it is also essential for TEB maintenance and mammary epithelial cell self-renewal (Lu et al., 2008; Parsa et al., 2008; Pond et al., 2013). The ligands involved are less well defined. Fgf10 is the main ligand expressed in the mammary stroma (Zhang et al., 2014), but assessing its role in vivo has been challenging, as the germline deletion results in the absence of most mammary buds and the mice do not survive beyond birth (Mailleux et al., 2002). Pathways inhibiting mammary gland development have also been identified, with Tgfβ1 being the major pathway (Silberstein and Daniel, 1987). Studies on postnatal mammary glands indicate that Tgfβ1 controls branching mainly by suppressing proliferation (Ewan et al., 2002). On the other hand, in vitro studies suggest that branches might initiate at sites with a local minimum concentration of Tgfβ1, implying a specific role in branch patterning (Nelson et al., 2006).
Much of our current understanding of lung, kidney and salivary gland branching morphogenesis stems from studies using established ex vivo organ cultures (Goodwin and Nelson, 2020; Hauser and Hoffman, 2015; Lang et al., 2021; Short and Smyth, 2016). Tools to study mammary ductal growth and branching have been limited because the commonly used system in mammary gland research, the 3D organoid culture, is typically devoid of stromal cells and, hence, the crucial epithelial-stromal crosstalk (Shamir and Ewald, 2014). Some recent organoid cultures do include stromal cells, yet their composition may not fully recapitulate that observed in vivo (Hume et al., 2018). Although numerous mouse studies have reported pubertal phenotypes, embryonic branching morphogenesis has remained a somewhat neglected area of research, possibly owing to challenges in working with small mammary rudiments and/or the early developmental arrest observed in mice deficient in some key pathways, such as Fgf10/Fgfr2b and Wnt (Andl et al., 2002; Mailleux et al., 2002; Spina and Cowin, 2021).
In this study, we aimed to decipher the principles of branching pattern formation in the mammary gland by combining 3D morphometrics with an ex vivo embryonic mammary gland culture method that preserves the epithelial-mesenchymal tissue interactions (Voutilainen et al., 2013). Our results show that the topologies of embryonic mammary glands are highly variable, suggesting a stochastic branching process, yet the average branching frequency remains remarkably constant throughout development and is therefore ‘stereotypical’. Branch point formation, however, was sensitive to exogenous manipulation of Fgf10 and Tgfβ pathways. Time-lapse imaging revealed that lateral branching and tip bifurcation jointly drive embryonic mammary gland morphogenesis, a finding supported by the 3D optical projection tomography data. Our data also imply that tissue-level rearrangements may contribute to ductal elongation.
RESULTS
Embryonic mammary gland branches stochastically but at a constant rate
To characterize the branching of embryonic mammary glands in 3D, fixed mammary glands were collected and scanned with optical projection tomography (OPT). Samples were dissected at embryonic day (E) 16.5 (when the branching morphogenesis commenced in most anterior glands), E17.5 and before birth at E18.5 from mice of mixed background. In total, 101 samples of mammary glands 1 to 5, ranging from unbranched mammary primordia (MG5 at E16.5) to glands with up to 25 ductal tips (MG1 at E18.5) were analyzed (Fig. 1A; Fig. S1A). As previously suggested (Veltmaat et al., 2003), the anterior mammary primordia were noticeably larger than the posterior primordia (Fig. 1A and Movie 1), and the patterns of ramifications differed from embryo to embryo, as well as between the left and right side of the same embryo (Fig. S1B). At E16.5, the 3D structure of the mammary gland was already relatively flat, reflecting its growth between the skin and muscle tissue (Fig. 1B, Movie 1).
Embryonic mammary glands develop asynchronously. (A) Representative OPT images of K14-Cre;mT/mG mammary glands (MG) 1 to 5 stained with anti-GFP antibody at E16.5, E17.5 and E18.5. (B) 180° rotational views of mammary gland 1 at E16.5, E17.5 and E18.5. (C) The terminology used throughout the article. Three types of branches, including the three branch types (main duct, segment and tip), are indicated. The main duct is the part between the nipple and the first branch point, the segment is an ‘internal branch’ between two branch points, and the tips are terminal branches. Each segment is assigned a generation number relative to the main duct. Scale bars: 150 µm.
Embryonic mammary glands develop asynchronously. (A) Representative OPT images of K14-Cre;mT/mG mammary glands (MG) 1 to 5 stained with anti-GFP antibody at E16.5, E17.5 and E18.5. (B) 180° rotational views of mammary gland 1 at E16.5, E17.5 and E18.5. (C) The terminology used throughout the article. Three types of branches, including the three branch types (main duct, segment and tip), are indicated. The main duct is the part between the nipple and the first branch point, the segment is an ‘internal branch’ between two branch points, and the tips are terminal branches. Each segment is assigned a generation number relative to the main duct. Scale bars: 150 µm.
For quantitative spatial analysis, images were analyzed with the Tree Surveyor software that enables derivation of network properties such as tip number, branch length and angle, tree volume, and branching generations (Short et al., 2013) (Fig. S1C). Quantification confirmed that branching increased from E16.5 to E18.5 based on both tip and generation numbers (Fig. 2A,B). In principle, the generation number can be used to infer the sequence of branching events, counting from the main duct of the ductal tree (for terminology used in the study, see Fig. 1C). When tip number was plotted against the number of generations per gland, it was evident that tips were not propagated, as expected in a stereotypical bifurcating system where each terminal branch gives rise to two new branches (Fig. 2C).
Embryonic mammary glands branch at a constant frequency. (A,B) Quantification of the number of ductal tips (A) and branch generations (B). (C) The theoretical number of tips (log2) per number of generations in the gland (black line) if branching occurs by repetitive bifurcations only and none of the terminal branches abort, and the experimentally observed relationship between the number of tips and the maximum number of generations found in the embryonic mammary glands (gray line). (D,E) Quantification of the total length (D) and total volume (E) of the ductal tree in mammary glands (MGs) 1-5. A two-tailed unpaired Student's t-test was used to assess statistical significance. (F) The ratio of total length to tip number is constant over developmental stages and across mammary glands 1-5. Statistical significance was tested using one-way ANOVA. Data in A,B,D,E are mean±s.d. (G-I) The correlations of the total ductal length and tip number of all mammary glands (G), the branch frequencies (total length/number of tips) and the number of tips (H), and the total volume of the ductal tree and the number of tips (I) were assessed using Pearson's correlation coefficient (r).
Embryonic mammary glands branch at a constant frequency. (A,B) Quantification of the number of ductal tips (A) and branch generations (B). (C) The theoretical number of tips (log2) per number of generations in the gland (black line) if branching occurs by repetitive bifurcations only and none of the terminal branches abort, and the experimentally observed relationship between the number of tips and the maximum number of generations found in the embryonic mammary glands (gray line). (D,E) Quantification of the total length (D) and total volume (E) of the ductal tree in mammary glands (MGs) 1-5. A two-tailed unpaired Student's t-test was used to assess statistical significance. (F) The ratio of total length to tip number is constant over developmental stages and across mammary glands 1-5. Statistical significance was tested using one-way ANOVA. Data in A,B,D,E are mean±s.d. (G-I) The correlations of the total ductal length and tip number of all mammary glands (G), the branch frequencies (total length/number of tips) and the number of tips (H), and the total volume of the ductal tree and the number of tips (I) were assessed using Pearson's correlation coefficient (r).
Quantification of the total length and volume also confirmed the increase in size in each gland between E16.5 and E18.5, again the anterior glands being larger than the posterior glands (Fig. 2D,E). However, the ratio of total ductal length and the number of ductal tips showed no overt difference across the glands (Fig. 2F), indicating that, despite the difference in size, the branch point frequency is constant and, hence, each gland branches essentially the same way. Strikingly, when all glands were pooled for analysis, the total length of the ductal tree and the number of ductal tips showed a near-perfect positive correlation (Pearson's r=0.968) (Fig. 2G). On average, there was one tip for every ∼260 µm±40 µm (mean±s.d.) of ductal length, a value that was constant over the developmental stage of the mammary glands (using tip number as a proxy for developmental stage), although less developed glands showed much higher variation compared with more developed glands (Fig. 2H). A strong positive correlation was also observed between the total volume and number of tips (Pearson's r=0.8629) (Fig. 2I). Pairwise comparison of E18.5 left-right mammary glands revealed prominent within-embryo variability, evidencing the stochastic nature of the branching process (Fig. S1D). Taken together, our OPT images show that each embryonic mammary gland displays a unique non-stereotypical branching pattern. Remarkably, however, they display a steady average increase in branch point frequency.
Both lateral branching and terminal bifurcations drive pre-pubertal mammary gland branching
As shown above (Fig. 2C), the topology of embryonic mammary glands did not confine to a stereotypical bifurcating system. This can be either because some terminal branches cease to grow, as previously suggested to take place during puberty (Scheele et al., 2017), or, alternatively, branches may form by side branching, which changes the generation dynamics, as a generation would be divided in two ‘retrospectively’. As it is not possible to distinguish between these two options based on fixed samples, we used an in-house developed organ culture system to perform time-lapse imaging of ex vivo cultured embryonic mammary glands. This culture set-up recapitulates growth in vivo relatively faithfully (Speroni et al., 2017; Voutilainen et al., 2013), although branching frequency was somewhat higher ex vivo (Fig. S2A). E13.5 mammary buds expressing epithelial GFP (K14-Cre; mT/mG) were dissected with the surrounding mesenchyme, and imaging was started on day 3 of culture, when branching had been initiated. Images were captured every 4 h for up to 4 days (Fig. 3A, Movie 2). Time-lapse imaging revealed that lateral branch formation was the main mode of branching, covering ∼75% of all branching events (Fig. 3B,C). Of these, 63% formed in the tips (i.e. in terminal branches; Fig. 3B, upper panel shows an example of a lateral branch forming in a tip) and 37% in segments (i.e. between two existing branch points). All glands monitored also had bifurcation events (24% of branching events). Terminal trifurcations were also observed, although rarely (<2% of all events). Analysis of the types of branching events by binned time intervals did not reveal any apparent difference in the timing of bifurcation and lateral branching events over the course of the 4-day imaging session (Fig. S2B).
Embryonic mammary gland branches via side branching and terminal bifurcations. (A) Time-lapse series of ex vivo cultured K14-Cre;mT/mG mammary glands. Mammary buds were dissected at E13.5, imaging was started after 3 days of culture and continued for up to 96 h. N, nipple. Arrows indicate the direction of growth. (B) An example of a time-lapse video used to assess and quantify the types of branching events: lateral branching (L), and bi- and tri-furcations (BFs). (C) Quantification of the types of branching events (n=149 branching events from 15 mammary glands). (D) Images depict how the branch point and tip closest to an emerging bifurcation (top) and lateral branch (below) were quantified. (E) Quantification of the distance to the closest branch point of nascent lateral branches and branches formed by bifurcation. For lateral branches, contra- (CO) and ipsilateral (IP) branch points were analyzed separately. (F) Quantification of the parental branch length at the time when a new lateral branch emerged between two ipsilateral branches (IP-IP), between two contralateral branches (CO-CO), between a terminal branch and an ipsilateral side branch (END-IP), and between a terminal branch and a contralateral side branch (END-CO). (G) Quantification of the distance to the closest branch tip of nascent lateral branches or branches formed by bifurcation. Data in E-G are mean±s.d. A two-tailed unpaired Student's t-test (E,G) or a Mann–Whitney test (F) was used to assess statistical significance. Scale bar: 100 µm (A).
Embryonic mammary gland branches via side branching and terminal bifurcations. (A) Time-lapse series of ex vivo cultured K14-Cre;mT/mG mammary glands. Mammary buds were dissected at E13.5, imaging was started after 3 days of culture and continued for up to 96 h. N, nipple. Arrows indicate the direction of growth. (B) An example of a time-lapse video used to assess and quantify the types of branching events: lateral branching (L), and bi- and tri-furcations (BFs). (C) Quantification of the types of branching events (n=149 branching events from 15 mammary glands). (D) Images depict how the branch point and tip closest to an emerging bifurcation (top) and lateral branch (below) were quantified. (E) Quantification of the distance to the closest branch point of nascent lateral branches and branches formed by bifurcation. For lateral branches, contra- (CO) and ipsilateral (IP) branch points were analyzed separately. (F) Quantification of the parental branch length at the time when a new lateral branch emerged between two ipsilateral branches (IP-IP), between two contralateral branches (CO-CO), between a terminal branch and an ipsilateral side branch (END-IP), and between a terminal branch and a contralateral side branch (END-CO). (G) Quantification of the distance to the closest branch tip of nascent lateral branches or branches formed by bifurcation. Data in E-G are mean±s.d. A two-tailed unpaired Student's t-test (E,G) or a Mann–Whitney test (F) was used to assess statistical significance. Scale bar: 100 µm (A).
As most explants had already started to branch when we initiated our imaging, we could not analyze the type of the very first branching event. To assess this indirectly, we turned to our OPT dataset and compared the length of the main duct in mammary sprouts that had not yet branched with those that had branched only once. Quantification revealed that the length of the main duct was 40% shorter and the volume 54% smaller in rudiments with two tips compared with unbranched ones (Fig. S2C,D), a finding that was poorly compatible with a tip-splitting event, because a significant shortening of an existing branch is difficult to envision occurring by any means other than formation of a new lateral branch. A pairwise comparison of the tip lengths in glands with two tips only showed that, although in some glands they were of equal length, others displayed a notable difference (Fig. S2E). In conclusion, our live-imaging data establishes that pre-pubertal branching morphogenesis uses a combination of lateral branching and terminal bifurcation events. Fixed sample data imply that the mode of the very first branching event is also stochastic.
Neighboring tips constrain the formation of new branches
Having identified both lateral and terminal branching events, we next analyzed in detail where they occur within the epithelial network and whether their locations follow any trend. For this purpose, we measured the distance of incipient lateral and bifurcating terminal branches to their closest neighboring tips and branch points (Fig. 3D). Tip splitting occurred at a minimum distance of 55 µm (average 103 µm) from the last branch point (Fig. 3E), indicating that terminal branches need to elongate to a certain length and/or grow to a certain size before bifurcating. In contrast, the occurrence of lateral branches was not as limited by their distance to existing branch points, although which side other branches were did make a difference (see Fig. S2F for a schematic representation of the different types of branch points). Branch points on the opposite sides (contralateral, CO) were tolerated at closer distances (minimum 9 µm; average 44 µm) than branch points on the same side (ipsilateral, IP) (minimum 15 µm; average 80 µm), possibly due to the fact that branching is constrained in the z-axis both in vivo and ex vivo (Fig. 3E). Accordingly, the length of the parental branch was dependent on the orientation of the closest branches. When a new branch emerged between two existing branches or a branch point and a tip, the length of the parental segment was significantly longer if the branch was an ipsilateral branch (IP-IP/END) compared with a contralateral branch (CO-CO/END) (Fig. 3F). Collectively, these data suggest that the proximity to other tips, rather than branch points themselves, might limit the location of a new lateral branch. Indeed, neighboring tips were at a minimum distance of 31 µm (mean 90 µm) from an emerging lateral branch; similar findings were made with bifurcations (Fig. 3G). Thus, neighboring tips, rather than branch points (except for ipsilateral branches), spatially constrain the emergence of new branches.
Spatial analysis of branch parameters predicts the pattern of branching, elongation and growth
Previous studies on the pubertal mammary gland localized branching and elongation to terminal branches, indicating that the length of a branch (segment) is determined once the tip bifurcates and is preserved thereafter (Scheele et al., 2017). To gain more insight into branch dynamics, we turned to our OPT dataset and separately analyzed terminal branches (tips), segments (branches between two branching points) and main ducts. Segments were found to be significantly shorter than terminal branches (Fig. 4A), a finding that cannot be explained if branching occurs via tip bifurcation only. However, as shown by our time-lapse data (Fig. 3B), both segments and tips could be shortened through lateral branching. The main ducts, on the other hand, were exceptionally long compared with both segments and terminal branches (Fig. 4A), as well as thick and large, as shown by the median diameter and volume, respectively (Fig. 4B,C; see also below). To further investigate how growth contributes to branch elongation in different parts of the gland, we visualized the dimensions of terminal branches and segments as a distribution of their length to the median diameter. The segments were narrower the longer they were, but the opposite was observed in terminal branches, suggesting that they may widen before they branch (Fig. 4D). We also assessed whether there was any correlation between the median diameter of the branch and the developmental stage of the gland (using tip number as a proxy). The segments, as well as the terminal branches, became narrower as the gland advanced in development (Fig. 4E).
Branch dynamics in embryonic mammary glands. Branch parameters were quantified with Tree Surveyor. (A-C) Quantification of the length of main ducts (n=89), segments (n=409) and terminal branches (n=592) (A), the median diameter of the length of main ducts (n=73), segments (n=344) and terminal branches (n=465) (B), and the volume of the main ducts (n=73), segments (n=347) and terminal branches (n=496) (C) from all mammary glands at all developmental stages. The data are median (line) with 25th and 75th percentiles (hinges) plus the min to max ranges (whiskers). Statistical significance was tested with one-way ANOVA with a Dunn's post-hoc comparison. (D,E) The correlation of the median diameter to branch length (D) and the number of tips (E) of segments (n=344) and terminal branches (n=465) was assessed using Pearson's correlation coefficient (r). (F) Length distributions of segments and terminal branches. Statistical significance was tested using a Kolmogorov–Smirnov test. (G) Correlation of the average segment length of each mammary gland to the number of tips (n=35) was assessed using Pearson's correlation coefficient (r). (H,I) Schematics depicting the local ‘bifurcation’ (n=491) (H) and dihedral (n=395) (I) branch angles, and their quantifications. Rayleigh test for nonuniformity, H0=random (for main duct parameters, only branched glands were included).
Branch dynamics in embryonic mammary glands. Branch parameters were quantified with Tree Surveyor. (A-C) Quantification of the length of main ducts (n=89), segments (n=409) and terminal branches (n=592) (A), the median diameter of the length of main ducts (n=73), segments (n=344) and terminal branches (n=465) (B), and the volume of the main ducts (n=73), segments (n=347) and terminal branches (n=496) (C) from all mammary glands at all developmental stages. The data are median (line) with 25th and 75th percentiles (hinges) plus the min to max ranges (whiskers). Statistical significance was tested with one-way ANOVA with a Dunn's post-hoc comparison. (D,E) The correlation of the median diameter to branch length (D) and the number of tips (E) of segments (n=344) and terminal branches (n=465) was assessed using Pearson's correlation coefficient (r). (F) Length distributions of segments and terminal branches. Statistical significance was tested using a Kolmogorov–Smirnov test. (G) Correlation of the average segment length of each mammary gland to the number of tips (n=35) was assessed using Pearson's correlation coefficient (r). (H,I) Schematics depicting the local ‘bifurcation’ (n=491) (H) and dihedral (n=395) (I) branch angles, and their quantifications. Rayleigh test for nonuniformity, H0=random (for main duct parameters, only branched glands were included).
Branch length measures elongation between two branching events that may or may not be sequential, as we have learned. In the pubertal mammary gland, the timing between consecutive branching events appears to be random (Hannezo et al., 2017). To understand how the spacing of branching events is established, we analyzed the distribution of segments and terminal branch lengths, reflecting the probability to branch. The distribution of segments was skewed towards shorter lengths, indicating that there is a tendency to minimize spacing between two branch points (Fig. 4A,F; Fig. S3A). A notable fraction (∼25%) of all segments was shorter than 55 µm, which is the minimum distance from a bifurcating tip to the closest branch point (Fig. 3D), a finding further supporting our conclusions on the contribution of side branching to mammary gland morphogenesis. Given that spatial constraints become greater the shorter the segment is, this minimum likely approaches the maximum epithelial packing density that is tolerated. Interestingly, when the distribution of the average segment length of each gland was analyzed against developmental time, we noticed much higher variation in the younger glands compared with the older ones (Fig. 4G).
The branching angles also contribute to the topology of the epithelial network. Bifurcating TEBs branch with an average angle of 70-75°, thereby biasing the global orientation of the ductal tree along the long axis of the fat pad (Nerger et al., 2021; Paine et al., 2016). We found the local ‘bifurcation’ angles of branch points (the angle between the start direction of two branches) to be larger, focusing around 90° (Fig. 4H). This finding is in line with the notion that the longitudinal bias is not observed in 2-week-old pre-pubertal mice (Nerger et al., 2021), In contrast, rotation between consecutive branch points (dihedral angle) was small (Fig. 4I), consistent with the rather two-dimensional architecture of the mammary gland.
Tgfβ1 and Fgf10 change branch point frequency
Mammary branching morphogenesis is regulated by conserved secreted signaling molecules, but their contribution to growth versus branch patterning is challenging to discern (Macias and Hinck, 2012; Myllymäki and Mikkola, 2019; Paine and Lewis, 2017). Some may simply modulate growth, whereas others might affect the propensity to branch. Our ex vivo culture system provides a means to distinguish between these two options, prompting us to assess how Tgfβ1 and Fgf10 regulate ductal morphogenesis. As there is no evidence of long-range guidance cues in the mammary gland stroma, and Fgf10 is uniformly expressed in the fat pad (Parsa et al., 2008), we reasoned that application of the candidate factors in the growth medium would best mimic the in vivo situation. Mammary buds were isolated from E13.5 embryos expressing epithelial GFP, allowed to grow for 2 days before applying either the signaling molecules or their inhibitors and then cultured for a further 5 days. We quantified the total ductal length and the tip number, and calculated their ratio to assess whether branching frequency was affected.
To assess the impact of Tgfβ1 on branching rate, we opted for a low concentration, given its well-known strong anti-proliferative effect (Moses et al., 2016; Silberstein and Daniel, 1987). When Tgfβ1 (2.5 ng/ml) was administered to explant cultures, a significant reduction in the number of tips was observed, and the total ductal length per tip number was increased, indicating a reduction in branching frequency (Fig. 5A; Fig. S4A). Treating explants with the TgfβRI inhibitor SB431542 had the opposite effect. The number of tips and total ductal length both increased, yet without a change in their ratio (Fig. 5B; Fig. S4B). To complement these experiments, we administered Tgfβ1 locally on E16.5 mammary explants in which some branches had already formed. Beads soaked with Tgfβ1, or BSA as a control, were added on day 1 (E16.5+1 day), explants were cultured for additional 2 days and the distance of the closest branch to the bead was measured at both time points. Although the explants with control and both Tgfβ1 beads grew and branched, the movement of the nearby branches towards the Tgfβ1 bead was significantly impaired (Fig. 5C,D).
Ex vivo manipulation of cultured mammary glands with Tgfβ1 and Fgf10 pathways. (A,B,E,F) Mammary buds were dissected at E13.5, the indicated signaling molecules/pathway inhibitors were added on day 2 and the explants were cultured for a further 5 days. The number of tips and total ductal length and their ratio were quantified on day 7. (A) Tgfβ1 (2.5 ng/ml) (nctrl=15, nTgfβ1=20 from two independent experiments). (B) SB431542 (10 µM) (nctrl=60, nSB431542=65 from eight independent experiments). (E) Fgf10 (100 ng/ml) (nctrl=21, nFgf10=23 from three independent experiments). (F) SU5402 (2.5 µM) (nctrl=24, nSU5402=21 from three independent experiments). (C,G) Mammary glands were dissected at E16.5 and beads (outlined in red) soaked with BSA (C,G), Tgfβ1 (C) or Fgf10 (G) were added on day 1 and the explants were cultured for additional 2 days. (D,H) Quantification of the relative distance (day1-day3/day1) between the closest tip and the bead for (D) Tgfβ1 (50 µg/ml) (nctrl=45, nTgfβ1=45 from four independent experiments) and (H) Fgf10 (100 µg/ml) (nctrl=38, nFgf10=38 from four independent experiments). Data are mean±s.d. Statistical significance was assessed using a two-tailed unpaired Student's t-test (B, total length/number of tips; D) or a Mann–Whitney test [A,B (number of tips and total ductal length),E,F,H]. Scale bars: 200 µm.
Ex vivo manipulation of cultured mammary glands with Tgfβ1 and Fgf10 pathways. (A,B,E,F) Mammary buds were dissected at E13.5, the indicated signaling molecules/pathway inhibitors were added on day 2 and the explants were cultured for a further 5 days. The number of tips and total ductal length and their ratio were quantified on day 7. (A) Tgfβ1 (2.5 ng/ml) (nctrl=15, nTgfβ1=20 from two independent experiments). (B) SB431542 (10 µM) (nctrl=60, nSB431542=65 from eight independent experiments). (E) Fgf10 (100 ng/ml) (nctrl=21, nFgf10=23 from three independent experiments). (F) SU5402 (2.5 µM) (nctrl=24, nSU5402=21 from three independent experiments). (C,G) Mammary glands were dissected at E16.5 and beads (outlined in red) soaked with BSA (C,G), Tgfβ1 (C) or Fgf10 (G) were added on day 1 and the explants were cultured for additional 2 days. (D,H) Quantification of the relative distance (day1-day3/day1) between the closest tip and the bead for (D) Tgfβ1 (50 µg/ml) (nctrl=45, nTgfβ1=45 from four independent experiments) and (H) Fgf10 (100 µg/ml) (nctrl=38, nFgf10=38 from four independent experiments). Data are mean±s.d. Statistical significance was assessed using a two-tailed unpaired Student's t-test (B, total length/number of tips; D) or a Mann–Whitney test [A,B (number of tips and total ductal length),E,F,H]. Scale bars: 200 µm.
The addition of Fgf10 (100 ng/ml) to mammary gland explant cultures significantly increased both the number of ductal tips and total ductal length (Fig. 5E; Fig. S4C), the latter indicating that Fgf10 enhances growth, as expected. In addition, branching frequency was augmented, as shown by the analysis of total ductal length per number of tips (Fig. 5E). Conversely, the broad-spectrum Fgf receptor inhibitor SU5402 severely hindered growth, as evidenced by the reduced total length of the ductal network and tip number, although branch frequency remained unaffected (Fig. 5G; Fig. S4D). We also applied Fgf10 locally but observed no difference in branch elongation toward the control and Fgf10-soaked beads (Fig. 5G,H).
To gain insights on the mechanisms by which Tgfβ1 and Fgf10 impact growth and/or branching, we assessed their effect on cell cycle dynamics in E16.5+3 day cultured mammary glands isolated from the constitutive bicistronic Fucci-2a cell cycle reporter mouse, where cells in the S/G2/M and G1/G0 phase express green and red nuclear fluorescent reporters, respectively (Mort et al., 2014; Lan et al., 2024). Supplementation of Tgfβ1 (2.5 ng/ml) significantly reduced the proportion of cells in the S/G2/M phase, whereas exogenous Fgf10 (100 ng/ml) had no apparent effect (Fig. S4E,F).
Tissue rearrangements contribute to main duct elongation
Next, we wanted to better understand the process of branch elongation and first focused on the main duct, which we can easily identify in fixed glands by its connection to the skin. Unlike terminal branches and segments, main duct lengths conformed to a rather normal distribution (Fig. S3A,B), suggesting that branching is prevented close to the skin. Embryonic stage-by-stage analysis revealed statistically significant elongation of the main duct (Fig. 6A), a finding confirmed by plotting the main duct length against the developmental stage of the mammary gland (using the total tip number as a proxy) (Fig. 6B). This result was surprising, given that branch elongation has mainly been associated with the growing tips (Scheele et al., 2017), and prompted us to examine cell proliferation in the main duct before and after the onset of branching morphogenesis. To this end, we analyzed fixed E16.5 mammary glands (ranging in size from non-branched ones to those with five tips) of the Fucci cell cycle indicator mice (Sakaue-Sawano et al., 2008) (Fig. 6C). Most cells in the main duct were quiescent, regardless of the developmental stage, whereas the tips were proliferative, as expected (Fig. 6D), suggesting that mechanisms other than proliferation drive main duct elongation.
The main duct elongates and narrows as development advances. (A) Quantification of main duct length from mammary glands 1-5 at E16.5 (n=18), E17.5 (n=36) and E18.5 (n=35). The data are median (line) with 25th and 75th percentiles (hinges) plus the min to max ranges (whiskers). (B) Correlation of the main duct length to the number of tips of all glands was assessed using Pearson's correlation coefficient (r). (C) A representative image of an E16.5 MG3 from Fucci mice. Cells in S/G2/M are in cyan; cells in G1/G0 are in purple. The dashed line indicates the epithelial-mesenchymal border identified by EpCAM immunostaining. (D) Quantification of the proportions of cells in G1/G0 and S/G2/M (or expressing both markers, in yellow) in main ducts (Md) and branches (Br) in mammary glands with one (unbranched) (n=2), two (n=3), three (n=4), or four or five tips (n=5). Statistical significance was assessed using a one-way ANOVA with Tukey's post-hoc comparison. Data are mean±s.d. (E) Quantification of main duct length from mammary glands 1-5 at E16.5, E17.5 and E18.5. The data are median (line) with 25th and 75th percentiles (hinges) plus the min to max ranges (whiskers). (F) Quantification of the median main duct diameter from mammary glands 1-5 at E16.5, E17.5 and E18.5 (number of main duct parameters, only branched glands were included). The data are median (line) with 25th and 75th percentiles (hinges) plus the min to max ranges (whiskers). (G) Average main duct diameter plotted against the number of tips (n=84) reveals that the diameter gradually settles at ∼70 µm. (H) Average gland diameter plotted against the number of tips (n=71) reveals that the diameter gradually settles at ∼70 µm. (I-K) Representative images of branches (last segments before the tips) (I) whose width (J) and length (K) (n=21) were quantified from the time-lapse data at timepoints 0 (T=0), 12 and 28. T0 indicates the time when the measurement starts. Data are normalized to the initial width and length of each branch. Data are mean±s.d. Scale bars: 100 µm.
The main duct elongates and narrows as development advances. (A) Quantification of main duct length from mammary glands 1-5 at E16.5 (n=18), E17.5 (n=36) and E18.5 (n=35). The data are median (line) with 25th and 75th percentiles (hinges) plus the min to max ranges (whiskers). (B) Correlation of the main duct length to the number of tips of all glands was assessed using Pearson's correlation coefficient (r). (C) A representative image of an E16.5 MG3 from Fucci mice. Cells in S/G2/M are in cyan; cells in G1/G0 are in purple. The dashed line indicates the epithelial-mesenchymal border identified by EpCAM immunostaining. (D) Quantification of the proportions of cells in G1/G0 and S/G2/M (or expressing both markers, in yellow) in main ducts (Md) and branches (Br) in mammary glands with one (unbranched) (n=2), two (n=3), three (n=4), or four or five tips (n=5). Statistical significance was assessed using a one-way ANOVA with Tukey's post-hoc comparison. Data are mean±s.d. (E) Quantification of main duct length from mammary glands 1-5 at E16.5, E17.5 and E18.5. The data are median (line) with 25th and 75th percentiles (hinges) plus the min to max ranges (whiskers). (F) Quantification of the median main duct diameter from mammary glands 1-5 at E16.5, E17.5 and E18.5 (number of main duct parameters, only branched glands were included). The data are median (line) with 25th and 75th percentiles (hinges) plus the min to max ranges (whiskers). (G) Average main duct diameter plotted against the number of tips (n=84) reveals that the diameter gradually settles at ∼70 µm. (H) Average gland diameter plotted against the number of tips (n=71) reveals that the diameter gradually settles at ∼70 µm. (I-K) Representative images of branches (last segments before the tips) (I) whose width (J) and length (K) (n=21) were quantified from the time-lapse data at timepoints 0 (T=0), 12 and 28. T0 indicates the time when the measurement starts. Data are normalized to the initial width and length of each branch. Data are mean±s.d. Scale bars: 100 µm.
Next, we analyzed the main duct volume and diameter at E16.5-E18.5. These data clearly show that although the main duct elongates, it does not grow in volume, but narrows in diameter (Fig. 6E-G). Together, these data show that the main duct elongates not by increasing its volume, but mainly by changing its aspect ratio from shorter and thicker to longer and thinner. The average diameter of not only the main duct but also the entire gland narrowed during branching morphogenesis, settling down at ∼70 µm (Fig. 6G,H; see also Fig. 4E), which is a value very close to that of mature ducts in pubertal glands (Paine et al., 2016). Unfortunately, for reasons that are not clear to us, the elongation of the main duct was rarely observed in our ex vivo cultures, precluding analysis by time-lapse imaging. Instead, we asked whether a similar process occurs elsewhere in the gland, and assessed whether segments could also narrow down and elongate even ‘behind’ a branching event (Fig. 6I). To this end, we measured the dimensions of newly formed segments from the time-lapse videos. Indeed, a significant change was observed in the aspect ratio: the segments elongated while their diameter narrowed (Fig. 6J,K).
Loss of Vangl2 compromises growth yet increases branch point frequency
Narrowing of the main duct without a change in the volume indicates convergent extension, a mechanism where the tissue converges (narrows) along one axis to lengthen it on a perpendicular axis, either via cellular movements or coordinated cell shape changes (Butler and Wallingford, 2017). Convergent extension is mediated through the conserved planar cell polarity (PCP) pathway, which in mammals is established by the asymmetrically distributed Van Gogh-like (Vangl) membrane proteins Vangl1 and Vangl2. To determine whether main duct elongation is regulated by PCP/Vangl2, we collected Vangl2-null and wild-type control C57Bl/6 mammary glands at E17.5, when branching had started in all glands and the main duct was elongating, and used 3D confocal microscopy for analysis (Fig. 7A; Fig. S5A). Again, anterior glands were larger than posterior glands, although in the C57Bl/6 background MG2 was the most developed gland (Fig. S5A,B). The phenotype of Vangl2-null mammary glands appeared variable (Fig. 7A; Fig. S5A), verified by the quantifications of tip number and total volume, and best exemplified by MG1-3 (Fig. S5B,C). For example, the number of tips was not affected except in MG1: Vangl2 mutants had significantly more tips, yet there was no change in the volume. On the other hand, MG2 and MG3 of Vangl2 mutants were greatly reduced in size although the number of tips was not altered (Fig. S5B,C). Together these findings suggest that Vangl2 may regulate the rate of branching. To test this hypothesis, we compared control and Vangl2 mutant mammary glands of equal tip numbers. The volume of Vangl2 mutant glands was consistently smaller (Fig. 7B) and hence the tip number per volume of the gland was significantly (∼60%) higher in Vangl2 mutants compared with their wild-type littermates (Fig. 7C), indicating higher propensity of branch point formation despite overall growth impairment. Analysis of the main ducts revealed that, on average, they were 35% smaller in volume, 26% shorter in length, but only 7% smaller in average diameter (Fig. 7D-F; Fig. S5D-F). Normalization of the main duct volume, length and average diameter to tip number confirmed that mainly the volume and length were affected, not the diameter (Fig. S5D-F).
Loss of Vangl2 compromises mammary gland growth but accelerates branching. (A) Representative confocal maximum intensity projection images of EpCAM-stained E17.5 mammary gland (MG) 2 of wild-type control (Vangl2+/+) and two Vangl2 knockout (KO) embryos. Scale bar: 100 µm. (B) Correlation of the mammary gland volume to number of tips in wild-type littermates and Vangl2 KO embryos was assessed using Pearson's correlation coefficient (r) (nwt=49 and nVangl2KO=57). (C) Quantification of the ratio of the number of ductal tips to the volume of glands of mammary glands 1-5 in wild-type littermates and Vangl2 KO embryos. Data are mean±s.d. (D-F) Main duct volume (D), main duct length (E) and main duct diameter (F) of mammary glands 1-5 in wild-type and Vangl2 KO embryos. Data are mean±s.d. (G) Representative confocal maximum intensity projection images of EpCAM-stained E18.5 mammary gland (MG) of wild-type control (Vangl2+/+) and Vangl2 knockout (KO) embryos. Scale bar: 50 µm. (H) Quantification of relative width of the main duct was measured near the first branch point/close to the nipple in wild-type and Vangl2 KO embryos (nwt=14 and nVangl2KO=15). Data are mean±s.d. (I) Cell segmentation based on EpCAM in the wild-type and Vangl2 KO mammary glands. Scale bar: 30 µm. (J) Quantification of the average cell volumes in wild-type (n=14) and Vangl2 KO (n=15) main ducts. Data are mean±s.d. (K) A schematic showing the cell orientation analysis in the main duct, and Rosa plots of cell orientation angles with respect to the elongating main duct; nwt=10,766 cells, nVangl2KO=22,008 cells. (L) Cell neighbor analysis in the E18.5 main duct of wild-type and Vangl2 knockout mammary glands (n=100 cells from 5 glands). Data are mean±s.d. Statistical significance was assessed using a two-tailed unpaired Student's t-test (H,J) or a Mann–Whitney test (C-F,L). Watson's U2 test and Rayleigh's test were used to test for nonuniformity, H0=random (K).
Loss of Vangl2 compromises mammary gland growth but accelerates branching. (A) Representative confocal maximum intensity projection images of EpCAM-stained E17.5 mammary gland (MG) 2 of wild-type control (Vangl2+/+) and two Vangl2 knockout (KO) embryos. Scale bar: 100 µm. (B) Correlation of the mammary gland volume to number of tips in wild-type littermates and Vangl2 KO embryos was assessed using Pearson's correlation coefficient (r) (nwt=49 and nVangl2KO=57). (C) Quantification of the ratio of the number of ductal tips to the volume of glands of mammary glands 1-5 in wild-type littermates and Vangl2 KO embryos. Data are mean±s.d. (D-F) Main duct volume (D), main duct length (E) and main duct diameter (F) of mammary glands 1-5 in wild-type and Vangl2 KO embryos. Data are mean±s.d. (G) Representative confocal maximum intensity projection images of EpCAM-stained E18.5 mammary gland (MG) of wild-type control (Vangl2+/+) and Vangl2 knockout (KO) embryos. Scale bar: 50 µm. (H) Quantification of relative width of the main duct was measured near the first branch point/close to the nipple in wild-type and Vangl2 KO embryos (nwt=14 and nVangl2KO=15). Data are mean±s.d. (I) Cell segmentation based on EpCAM in the wild-type and Vangl2 KO mammary glands. Scale bar: 30 µm. (J) Quantification of the average cell volumes in wild-type (n=14) and Vangl2 KO (n=15) main ducts. Data are mean±s.d. (K) A schematic showing the cell orientation analysis in the main duct, and Rosa plots of cell orientation angles with respect to the elongating main duct; nwt=10,766 cells, nVangl2KO=22,008 cells. (L) Cell neighbor analysis in the E18.5 main duct of wild-type and Vangl2 knockout mammary glands (n=100 cells from 5 glands). Data are mean±s.d. Statistical significance was assessed using a two-tailed unpaired Student's t-test (H,J) or a Mann–Whitney test (C-F,L). Watson's U2 test and Rayleigh's test were used to test for nonuniformity, H0=random (K).
To gain deeper insights into the role of Vangl2 in main duct elongation, we focused on embryonic stage E18.5, the latest stage we could obtain live mutant embryos. Mammary glands were stained with EpCAM and phalloidin, and the main ducts were imaged at high resolution in 3D. Again, the main ducts of Vangl2 mutants were abnormally short and, in contrast to wild type, often ‘V-shaped’, i.e. narrower near the nipple and wider near the first branch point (Fig. 7G), a finding confirmed by quantification (Fig. 7H). As two distinct processes take place during late embryogenesis, lumen formation and ductal elongation, and might affect one another, we first analyzed the former. In both control and Vangl2 mutants, formation of microlumini was evident and associated with high levels of F-actin, suggesting no major impairment in lumen formation (Fig. S5G). For cell level analysis, cells were segmented in 3D based on EpCAM staining (Fig. 7I). There was no difference in average cell volume (Fig. 7J), but there was a significant (although minor) difference in cell sphericity (Fig. S5H) between the two genotypes.
Convergent extension can be either passive (the surrounding environment deforms the tissue) or active (cells of the elongating tissue generate the necessary forces) (Belmonte et al., 2016). Cell intercalation takes place in both, but typically the axis of cell elongation is parallel to the axis of tissue elongation in the former and perpendicular in the latter. Quantification of the longest axis of the cell with respect to the axis of the elongating main duct revealed that cells in wild-type main duct were non-randomly distributed, and preferentially oriented perpendicular to the long axis of the main duct (Fig. 7K: 0° equals to perfect alignment with the elongation axis). The same was true in Vangl2 mutants. However, there was a statistically significant difference in the cell orientation between the two genotypes; however, the difference was very small and therefore potentially not biologically relevant. Cell neighbor analysis in the main duct revealed that Vangl2 KO cells have significantly fewer neighboring cells compared with the cells in wild-type ducts (Fig. 7L), which is indicative of compromised cell-cell interactions that may contribute to the observed phenotype.
DISCUSSION
In this study, we took advantage of the Tree Surveyor software designed to interrogate branching morphogenesis in 3D (Short et al., 2013) and ex vivo time-lapse imaging (Voutilainen et al., 2013) to provide a comprehensive description of the embryonic mammary gland branching morphogenesis by mapping and quantifying the embryonic ductal network and its dynamic changes over time. The topologies of the individual epithelial trees were highly variable, as reported for pubertal and adult glands (Hadsell et al., 2015; Scheele et al., 2017), evidencing that mammary branching is stochastic from the very beginning. Our findings are in line with the idea of mammary branching morphogenesis as ‘random branching and exploratory walk’, ultimately leading to efficient filling of the available space (Hannezo et al., 2017; Nerger et al., 2021). In pubertal glands, the model entails termination of branches due to crowding, i.e. when tips arrive too close to another duct, tip or the border of the fat pad (Hannezo et al., 2017). However, such annihilation behavior is likely to be less relevant (or non-existing) in the embryo, where the stroma can readily accommodate much larger ductal trees, as exemplified by the K14-Eda mouse model with a fivefold increase in the number of ductal tips (Voutilainen et al., 2012). In the embryonic gland, branch angles were relatively homogenous, centering around 90°, suggesting a means to maximally avoid the parental duct (in the case of lateral branching) or the sister tip (after bifurcation), thereby leading to an efficient ‘spread’ of the epithelial network in a yet unobstructed stroma. The non-stereotypy of the ductal trees was most evident in the diverse branch lengths; however, these did not follow a random distribution but gravitated toward the shorter end with a subset of very long branches, thereby generating highly heterogeneous epithelial networks. Despite these polymorphic branch patterns, branching was, on average, remarkably constant, as shown by the invariant ratio of the total length of the epithelial tree to branch number, and hence in this respect remarkably ‘stereotypical’.
Molecular regulation of branching morphogenesis has been investigated extensively (Paine and Lewis, 2017; Spina and Cowin, 2021). Many mouse models display impaired ductal morphogenesis, yet in many cases it remains unclear whether growth (the size of the ductal tree) and/or branching (the branch point frequency per ductal length) is affected, as few studies report these parameters rigorously. In addition, stroma-free organoid assays have provided valuable insights (Goodwin and Nelson, 2021; Huebner et al., 2016; Nelson et al., 2006), although to what extent these findings reflect in vivo behavior of the mammary epithelium is uncertain. The autocrine growth-inhibitory function of Tgfβ1 on postnatal mammary morphogenesis is well described (Ewan et al., 2002; Nelson et al., 2006; Silberstein and Daniel, 1987), but its role during embryogenesis has remained elusive. Our ex vivo organ culture experiment showed that ductal tips avoided locally released Tgfβ1 and that culture medium supplemented with Tgfβ1 decreased branching frequency. This occurred at a low concentration of Tgfβ1 that did not impair overall growth, suggesting a specific effect on branching probability. In contrast, suppression of Tgfβ activity ex vivo promoted growth, but did not alter branching frequency. These data are in line with previous analysis of pubertal Tgfb1+/- mice expressing ∼10% of wild-type levels of Tgfβ1, which show accelerated ductal growth without any apparent change in branch patterning (Ewan et al., 2002). Collectively, these data imply that endogenous Tgfβ1 signaling mainly functions to limit growth. Yet there may be a narrow physiological range of Tgfβ1 levels and/or local variation in its activity that contribute to the choice of branch initiation, as suggested by in vitro studies (Pavlovich et al., 2011).
Multiple in vivo studies link the Fgf pathway to mammary morphogenesis. In a genetic mosaic model, Fgfr2-null epithelial cells are outcompeted by the wild-type cells in TEBs due to a proliferation defect, and conditional attenuation of Fgfr1/2 activity severely compromises ductal growth (Lu et al., 2008; Parsa et al., 2008; Pond et al., 2013), establishing an important function in proliferation, but revealing less about contribution to branching. The role of Fgfr ligands is not well understood either. In 3D organoid culture, Fgf10 fails to induce branching, unlike Fgf2, which is the commonly used ‘branching inducer’ in these assays, but does function as a chemoattractant (Zhang et al., 2014); however, we found no evidence for such a role in the context of the whole explant. Instead, when added to the culture medium, exogenous Fgf10 enhanced not only growth, but also branching. However, we observed no change in cell cycle progression, implying that the effect of Fgf10 is mediated by other cellular behaviors, such as survival and/or motility. On the other hand, attenuated Fgf activity decreased ductal growth without impairing branching frequency. One possible explanation for this counterintuitive finding is that some other signaling pathway may compensate for the initiation of branch formation (but not growth, as shown in vivo) when Fgf activity is suppressed. Alternatively, proliferation and branch point formation may require different Fgf activity levels and/or even distinct downstream pathways.
Studies from other branching organs suggest that tube elongation may be driven via multiple mechanisms, including convergent extension, collective migration and oriented cell divisions (Kunimoto et al., 2017; Tang et al., 2018; Torban and Sokol, 2021; Varner and Nelson, 2014). The PCP pathway regulates all these processes and has been implicated in tubule elongation in the kidney, where mutations in the PCP components, including Vangl2, result in dilated ducts (Kunimoto et al., 2017; Torban and Sokol, 2021). Growth of the mammary ductal network is undoubtedly fueled by the proliferative tips (Lu et al., 2008), and we have recently identified oriented cell motility as an important cellular mechanism driving ductal growth (Myllymäki et al., 2023). Our current data suggest that other mechanisms may also contribute to ductal elongation: a mere change in the ductal aspect ratio can lead to elongation, implying involvement of tissue rearrangements. Cells of the main duct were oriented perpendicular to the elongating main duct, suggesting that this might be an active process; however, testing this hypothesis must await the development of a culture set-up capable of capturing this phenomenon ex vivo. A recent study assessing the mammary phenotype of mice homozygous for the Looptail (Lp) allele encoding a dominant-negative Vangl2 protein that also interferes with the related Vangl1 protein (Yin et al., 2012) reported that 30% of Lp mammary epithelia transplanted into cleared fat pads generated outgrowths with wide ducts (Smith et al., 2019). This phenotype was not observed in single conditional Vangl1 or Vangl2 mutants, but was shared with mice heterozygous for Prickle2, another PCP component (Smith et al., 2019), indicating that PCP is also involved in ductal elongation in the mammary gland. Our data show that Vangl2-null embryos had shorter main ducts, but without a concurring increase in diameter, as one might expect if convergent extension was compromised. Instead, branching frequency was increased. This may indicate that the cell movements/rearrangements necessary for elongation are only partially impaired, e.g. due to redundancy with the related Vangl1 (Smith et al., 2019). Alternatively, some other Vangl2-dependent processes, such as directional migration, might be impaired. Further studies are warranted to clarify the role of Vangl2/PCP in mammary ductal morphogenesis.
Thus far, lateral branching has been reported to occur only in developing lungs and mammary glands. In the latter, side branches emerge during the estrous cycle and early pregnancy under the influence of hormonal cues, mainly progesterone (Brisken and Ataca, 2015). Here, by using time-lapse imaging we unequivocally establish that side branching is the predominant mode of branching in the embryonic mammary gland, although tip bifurcations are also common. Our data indicate that the proximity to neighboring tips, rather than branch points, spatially constrain the emergence of new lateral branches. Whether lateral branching takes place during puberty has been a matter of controversy. Pubertal glands are characterized by long extending ‘primary’ ducts that repeatedly split into two new ducts by bifurcation. Additionally, shorter ‘secondary’ branches stem from the primary ducts. Earlier research has commonly interpreted (many of) the shorter branches to result from side-branching events, a conclusion supported by whole mounts and histological sections, revealing the emergence of nascent branches from an existing duct, often close to its invading tip (Fata et al., 2004; Silberstein, 2001; Silberstein and Daniel, 1982; Sternlicht et al., 2006). A more recent study concluded that pubertal branching occurs almost exclusively through tip bifurcation (Scheele et al., 2017). Recent advances in intravital time-lapse imaging of postnatal mammary glands (Lloyd-Lewis, 2020) should facilitate the efforts to decipher whether lateral branching also occurs during the pubertal period or whether puberty is an exception in the journey of the developing mammary gland.
In conclusion, we report here the use of quantitative methods producing baseline data on embryonic mammary gland branching morphogenesis. The epithelial tree forms by a combination of side-branching and tip-bifurcation events (and to a small extent via trifurcations), similar to the developing lung (Metzger et al., 2008). In contrast to the lung, however, rotations between branching events are minimal, leading to the formation of a planar epithelial network from early on. In embryonic mammary glands, branch points form seemingly stochastically, leading to highly variable network topologies; yet development progresses predictably, owing to constant average branching frequency. How growth is ‘measured’ even when branching is stochastic and how, at the cellular level, different signaling pathways modulate the branching frequency remain fascinating challenges for future research.
MATERIALS AND METHODS
Mice
The following mouse strains were used: ROSAmT/mG (Jackson Laboratory stock no. 007576) (Muzumdar et al., 2007) (hereafter mT/mG mice) in a mixed C129SvJ/ICR background, K14-Cre43 (hereafter K14-Cre) (Andl et al., 2004) and Fucci (RIKEN Tsukuba BioResource Center) (Sakaue-Sawano et al., 2008) mice in an outbred NMRI background (Janvier Labs), and constitutive Fucci-2a mice in C57Bl/6 background (Mort et al., 2014; Lan et al., 2024). A Vangl2-null allele in C57Bl/6 background was derived by serially crossing Vangl2 ex2-3 floxed mice (Jackson Laboratory stock no. 025174) (Copley et al., 2013) with PGK-Cre mice universally expressing Cre (Jackson Laboratory stock no. 020811) (Lallemand et al., 1998) to permanently delete exons 2 and 3. All Vangl2 KO embryos displayed severe neural tube defects, confirming successful deletion (Copley et al., 2013). The Vangl2-null allele was genotyped with 5′-CATTTGTCTGTTGCTGGGTAA-3′ and 5′-TGGACTGTGTCTTGCCTACTG-3′ primers. The appearance of the vaginal plug was considered as embryonic day (E) 0.5. The embryonic stage was further verified by limb morphogenesis and other external criteria.
Ethics statement
All mouse experiments were approved by the local ethics committee and the National Animal Experiment Board of Finland (licenses KEK16-021, KEK19-019, KEK22-014 and ESAVI/2363/04.10.07/2017). Mice were euthanized with CO2 followed by cervical dislocation.
Optical projection tomography
E16.5, E17.5 and E18.5 K14-Cre;mT/mG female embryos were decapitated and bisected, and internal organs removed before fixation (4% PFA in PBS overnight at +4°C). Individual mammary glands were dissected with a small amount of surrounding skin under the fluorescent microscope and immunohistochemically stained as previously described (Alanentalo et al., 2007). Chicken IgY anti-GFP (GFP-1010, 1:500, Aves Labs) was used as the primary antibody, and Cy™3 AffiniPure Goat Anti-Chicken IgY (103-165-155, Jackson ImmunoResearch Europe) was used as the secondary antibody. Samples were mounted in SeaPlaque agarose (50100, Lonza), dehydrated with methanol, cleared with BABB solution as described previously (Alanentalo et al., 2007), and scanned with Bioptonics OPT Scanner 3001 M with 6.7 µm camera pixel size and rotation step of 0.900. Scanned images were reconstructed with NRecon (Bruker). Mammary gland epithelium was roughly segmented by ImageJ (TrakEM package) and three-dimensional analysis was carried out using Tree Surveyor software (Short et al., 2013). Maximum diameter was determined for each sample and skeletonization errors were manually corrected. In Tree Surveyor software, all branching points were assigned as bifurcations. If two branches seemingly originated from the same branching point (segment lengths <7 µm), they were assigned as a trifurcation to avoid misconfiguration of branching generations. Terminal branches (tips) with a median diameter value of 0 µm were considered artifacts and excluded from the analysis. The average diameter of the gland was calculated by formula 2 × √([total volume]/(π × [total length])). The average diameter of the main duct was calculated with the same formula with [volume of main duct] and [length of main duct]. For branch point angles, the local branch angle values were reported.
Immunofluorescence and confocal microscopy
For whole-mount immunofluorescence staining, the ventral skin containing mammary glands was dissected from female E16.5, E17.5 or E18.5 mouse embryos. The samples were fixed in 4% PFA at 4°C overnight and then washed with PBS for 3-4 h. After permeabilization with 0.3% PBST (0.3% Triton X-100 in PBS) for 1-2 h at room temperature, E16.5 and E17.5 samples were blocked with blocking buffer (5% normal donkey or goat serum, 0.5% BSA and 0.5% Triton X-100 in PBS) for 1 h. E18.5 samples were blocked without the permeabilization step with blocking buffer (5% goat serum, 1% BSA and 0.3% TritonX-100 in PBS) for 1 h. All samples were incubated at 4°C with rat anti-mouse CD326 (EpCAM, BD Pharmingen, 552370, 1:1000) and 10 μg/ml Hoechst 33342 (Molecular Probes/Invitrogen) in a blocking buffer for 2-3 days. After washing with 0.3% PBST for 3-4 h, EpCAM was detected with an Alexa Fluor 488- or Alexa Fluor 647-conjugated secondary antibody (1:500, Molecular Probes/Invitrogen) at 4°C for 2 days. For E18.5 samples, Phalloidin-Alexa 568 was included (1:300, Invitrogen) to visualize filamentous actin. After washing in 0.3% PBST for 2-3 h, individual mammary glands with limited surrounding mesenchyme were dissected in PBS under a stereomicroscope and mounted within Vectashield (Vector Laboratories). The images were acquired in 1.10 µm (Fucci samples) or 1.40 µm (Vangl2-null and wild-type littermate samples) intervals by using a Zeiss LSM700 laser scanning confocal microscope with Plan-Apochromat 25×/0.8 immersion objective (Carl Zeiss Microscopy). E18.5 Vangl2 KO and their littermate control samples were cleared in a solution containing sucrose and glycerol overnight, and mounted, then images of main ducts were acquired with Leica Sp8 upright confocal microscope at 1.04 µm intervals with HC PL APO 20×/0.75 IMM CORR CS2 objective with glycerol. Images of growth factor-treated Fucci-2a explants were acquired with the same microscope and objective at 1.0 µm intervals.
Organ culture and time-lapse imaging
Mammary buds were dissected from E13.5 embryos and cultured in a Trowell-type culture setting at the liquid-air interface, as previously described (Lan et al., 2022; Voutilainen et al., 2013). All embryos of the preferred genotype were used because both females and males develop similarly until E14.0. After dissecting the ventrolateral skin encompassing the mammary area, explants were treated for 25-35 min for Dispase II (Roche, 1.25 U/ml) followed by a 4-8 min treatment with 0.75% pancreatin and 0.25% trypsin (Sigma-Aldrich) at room temperature. After a 30-40 min incubation on ice in DMEM supplemented with 10% (vol/vol) FCS, 2 mM L-glutamine and 20 U/ml penicillin-streptomycin (Gibco/Invitrogen), excess epithelium surrounding the mammary buds was microsurgically peeled off, resulting in explants with mammary primordia attached to an ample amount of mesenchyme. Mammary buds 2 and/or 3 were used for quantification in growth factor experiments. The culture medium consisted of a 1:1 mixture of DMEM with 2 mM l-glutamine and F12 (Ham's Nutrient Mix) (both from Life Technologies) and was supplemented with 10% (vol/vol) FCS (HyClone/Thermo Scientific), 20 U/ml penicillin-streptomycin (Gibco/Invitrogen) and ascorbic acid (100 mg/l). Explants were cultured at 37°C with 5% CO2.
To test the effect of growth factors and signaling pathway inhibitors, mammary glands from wild-type E13.5 NMRI embryos were cultured as described above. The sex of the embryos used in the experiments was not controlled. The following proteins and inhibitors were reconstituted according to the manufacturer's recommendations and used at the indicated final concentrations: Fgf10 (R&D, 345-FG, 100 ng/ml with 2 µg/ml heparin), SU5402 (Calbiochem, 572630, 2.5 µM), Tgfβ1 (R&D, 7666-MB, 2.5 ng/ml) and SB431542 (Tocris Bioscience, 301836-41-9 K, 10 µM). Cultured mammary glands begin branching on day 3, reflecting the in vivo onset of branching (Voutilainen et al., 2013). To specifically assess the effect of the signaling molecules/inhibitors on branching rather than initial outgrowth, vehicle (untreated control; 2 µg/ml heparin in 0.1% BSA for Fgf10 and 4 mM HCl in 0.1% BSA in PBS for Tgfβ1) or proteins/inhibitors were added on day 2 of culture; medium with the vehicle, protein or inhibitor was changed every other day. From the second day of culture, explants were imaged once per day with a Zeiss Lumar microscope equipped with an Apolumar S 1.2× objective. For the bead experiment, mammary explants were collected from E16.5 embryos by micro-dissecting mammary glands 1-3 with the surrounding mesenchyme and cultured in the same Trowell type set up for 3 days (Lan et al., 2022). Beads soaked with Fgf10 (100 µg/ml) or Tgfβ1 (50 µg/ml) and control beads soaked with BSA were added on the next day at a distance of ∼250 µm from the nearest tip. After adding the beads, explants were imaged once per day every day with a Zeiss Axio Zoom.V16 microscope at 50× magnification. To assess the effect of the growth factors on cell cycle progression during branching morphogenesis, mammary glands 1-3 were dissected from E16.5 female Fucci-2a embryos as above, and cultured for 3 days in media supplemented with the vehicle control, Tgfβ1 (2.5 ng/ml) or Fgf10 (100 ng/ml). Growth factors were supplemented on day 0, day 2 and day 3, 5 h before fixation of the samples.
Multi-position, automated time-lapse imaging of K14-Cre;mT/mG explants was used to assess the types of branching events (lateral branching, bi/trifurcations). Explants were dissected as described above, except that each explant had only mammary gland 3 to avoid any potential interference from the closely positioned mammary gland 2. Explants were placed on Transwell inserts (Costar, Corning) and cultured on six-well plates, allowing multiposition imaging. From day 3 to day 7 of culture, ductal trees were imaged with 3I Marianas wide-field microscope equipped with a 10×/0.30 EC Plan-Neofluar Ph1 WD=5.2 M27 at 37°C with 6% CO2. The medium was changed on day 2, at the onset of imaging (day 3), and once during imaging, on day 5. Images were acquired with a LED light source (CoolLED pE2 with 490 nm/550 nm) every 4 h.
Image analysis
K14-Cre;mT/mG mammary gland explants
The types of branching events were assessed by careful examination of subsequent time point images. Because of two-dimensional imaging, not all parts of the glands were in focus at all time points and occasionally the type of branching remained uncertain. Quantification of branching types was carried out only in glands in which at least 60% of branching events could be reliably determined. Time-lapse images were analyzed with ImageJ. Time series were registered with StackReg plug-in (Thevenaz et al., 1998), background subtracted (rolling ball r=200) and Gaussian (sigma=2) filtered and binarized with the local threshold (Li method). Segmentation errors were manually corrected (i.e. gaps in mammary epithelium and separation of falsely connected branches), and the fill holes function and Shape Smoothing plug-in were applied. Segmented images were skeletonized and skeletonization errors were manually corrected. Analyze Skeleton in 2D time series script (https://imagej.net/Analyze_Skeleton_2D_time_series) was used for measuring the skeletons. Epithelial area was measured from segmented images by ROI Manager in ImageJ and average ductal diameter was derived using the formula [total area]/[total length] for each time point.
Ductal diameter before bifurcation was measured by fitting an ellipse perpendicular to the spline of the duct in the ductal end in the time point preceding the bifurcation event and measuring the major axis of the ellipse. Distances to the closest tip (bifurcations and lateral branching) and closest ipsilateral tip (lateral branching only) were measured using the ‘Straight line’ tool and ROI Manager in ImageJ. Distances to neighboring branch points and length of the parental segment were measured from the skeletonized images in ImageJ. The minimum ductal diameter of the parental segment after the branching event was approximated by fitting a ‘Rotated Rectangle’ to the narrowest point in the segmented images. The diameter was measured from each time point as three independent replicates that were inside a 2% difference. Measurements were not performed on branches forming next to each other at the same time point.
Ductal elongation measurement was carried out with ImageJ (Fiji) using the plug-in Analyze-skeleton-Analyze skeleton 2D/3D. The Excel file with branch length information was saved, x and y positions of the fragment of interest from the tagged skeleton image were found, and the branch length values corresponding to the x and y values were noted.
Effect of growth factors and/or inhibitors on ductal morphogenesis in ex vivo cultured mammary glands
Quantification of the number of branches, total ductal length of the mammary gland and, in the bead experiments, the distance of the nearest tip to the bead (which was placed at a distance of ∼250 µm from the nearest tip) were carried out with ImageJ (Fiji). The distance to the bead was measured on day 1 (after applying the bead) and on day 3, and data are presented as a relative distance of day 3 to day 1/day 1 (i.e. a value of 0 indicates no growth towards the bead; a value of 1 indicates that the tip had reached the bead).
Cell cycle status quantification
To quantify cell cycle dynamics of epithelial cells in the developing mammary gland, E16.5 female embryos expressing the Fucci transgenes were stained with EpCAM as described above to discriminate the mammary epithelium and mesenchyme. The acquired confocal images were analyzed manually with Zen software (V2.3, Carl Zeiss Microscopy). Briefly, optical sections were selected every 10-11 µm in the z-axis, and mKO2 (cells in G1/G0) or mAG (cells in S/G2/M) expressing mammary epithelial cells (EpCAM+) were counted manually in selected optical sections. The main duct area was manually separated in Imaris at the position where the first branching occurs, and cells were counted accordingly. To assess the effect of Tgfβ1 and Fgf10 on cell cycle status in E16.5+3d cultured mammary glands, female embryos of the bicistronic Fucci-2a mice were used instead, as the transgenes of the Fucci mice were silenced at late embryonic stages. For cell proliferation analysis, the 3D surface rendering was carried out based on EpCAM staining, which was used to mask the mesenchymal Fucci-2a reporters, allowing better visualization of the epithelial Fucci-2a signals. For quantification of Fucci-2a reporter-expressing cells, only tips (but not segments or the main duct) were included in the analysis, which was performed using the Spots function in Imaris software with 7 and 8 µm spot sizes for G1/G0 and M/G2/S cells, respectively.
Vangl2 knockout analysis
EpCAM whole-mount immunostained confocal images were processed with Imaris software for the 3D rendering of the mammary epithelium to visualize the ductal tree. The number of tips was counted manually based on the 3D view. Main duct analysis was carried out with the combined use of both Fiji ImageJ and Imaris software. The main duct area was manually separated, as described above. Masked stacks were made based on EpCAM staining and skeletonized with Fiji, and the main duct length was quantified with the plug-in script from Arganda-Carreras et al. (2010). Main duct volume was quantified from 3D rendered images in Imaris. The average main duct diameter was calculated using the formula 2 × √([total volume]/(π × [total length])).
For cell shape/volume analysis, the confocal image stacks were preprocessed with the ImageJ plug-in Noise2Void (Broaddus et al., 2020) with the N2V train and predict module to remove the background. The training was performed with 100 epochs, 200 steps per epoch, a batch size per step of 64 and a patch shape of 64. The neighborhood radius was adjusted to either 4 or 5, based on the quality of the images. Cells were then segmented in 3D using a cell membrane marker, EpCAM, with Imaris software version 10 (Bitplane). Segmentation was manually examined and poor-quality cells were removed from later analysis. Cell volume was obtained using Imaris. To assess cell orientation, a Reference Frame was manually adjusted to align with the long axis of the main duct for each sample. The vector presents the cell elongation direction (ellipsoid axis length C) and the coordinates of the cell center were extracted for analysis in R with package ‘circular’ (Agostinelli and Lund, 2017) and ‘tidyverse’ (Wickham et al., 2019). To measure the width of the main duct near the nipple and branch, we used the spot module of the Imaris software version 10 (Bitplane). Two spots were placed in 3D and the width of the duct was derived from the distance between the two spots and exported from Imaris.
Cell neighbor analysis
For cell neighbor analysis, 20 cells were randomly selected in two or three optical sections from the middle of the E18.5 main duct from five wild-type and five Vangl2 KO embryos, and the cell neighbors in contact with each cell were counted manually with ImageJ (Fiji).
Statistical analysis
All data were analyzed using Prism 7 and 9 (GraphPad Software), or by circular package for R (Agostinelli and Lund, 2017). The normality of the data was assessed with a Shapiro-Wilk test. Statistical tests used are indicated in the figure legends. P<0.05 was considered significant. Throughout the figure legends: *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001. Rayleigh's Z-test was used to assess whether cells of the main duct are randomly oriented with respect to the main duct by quantifying the angle between the longest aspect of the cell with respect to the elongating main duct (0°=parallel to the main duct; 90°=perpendicular to the main duct). Watson's U2 test was used to analyze the difference in angle distribution between the wild-type and Vangl2 KO samples.
Acknowledgements
We thank Ms Raija Savolainen, Ms Riikka Santalahti and Ms Agnes Viherä for excellent technical assistance, and all Mikkola lab members for stimulating discussions. Drs Ian Smyth and Kieran Short are acknowledged for sharing the Tree Surveyor software, and Dr Jukka Jernvall for advice. The mouse studies were carried out with the support of the HiLIFE Laboratory Animal Centre Core Facility, University of Helsinki. Confocal and wide-field microscopy and image analysis with Imaris were performed at the Light Microscopy Unit, Institute of Biotechnology, University of Helsinki, supported by HiLIFE and Biocenter Finland. Optical projection tomography was carried out at the BioImaging unit of the University of Oulu.
Footnotes
Author contributions
Conceptualization: M.L.M., J.P.S., R.L.; Methodology: J.P.S., R.L., S.-M.M., Q.L., R.P.-H., M.V., S.K., S.J.V.; Formal analysis: J.P.S., R.L., S.-M.M., Q.L.; Investigation: J.P.S., R.L., S.-M.M., Q.L., E.T., B.K.; Resources: M.L.M., S.J.V.; Writing - original draft: M.L.M., J.P.S., R.L., S.-M.M.; Writing - review & editing: M.L.M., J.P.S., R.L., S.-M.M., Q.L., S.K., S.J.V.; Visualization: J.P.S., R.L., S.-M.M., Q.L., E.T.; Supervision: M.L.M.; Project administration: M.L.M.; Funding acquisition: M.L.M.
Funding
This study was financially supported by the Research Council of Finland Center of Excellence Program (grants 272280 and 307421 to M.L.M.; and project grant 318287 to M.L.M.), by the Sigrid Juséliuksen Säätiö (M.L.M.), by the Jane ja Aatos Erkon Säätiö (M.L.M.), by the Syöpäsäätiö (M.L.M.) and by the Helsinki Institute of Life Science Fellow Program (M.L.M.). J.P.S. acknowledges support from the Suomen Kulttuurirahasto and from the Ella ja Georg Ehrnroothin Säätiö. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Data availability
All relevant data can be found within the article and its supplementary information.
Peer review history
The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.202179.reviewer-comments.pdf
References
Competing interests
The authors declare no competing or financial interests.