Oriented clonal cell dynamics enables accurate growth and shaping of vertebrate cartilage

*For correspondence: kaj.fried@ ki.se (KF); [email protected] (IA) Competing interests: The authors declare that no competing interests exist. Funding: See page 28 Received: 09 February 2017 Accepted: 16 April 2017 Published: 17 April 2017 Reviewing editor: Marianne Bronner, California Institute of Technology, United States This is an open-access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication. Oriented clonal cell dynamics enables accurate growth and shaping of vertebrate cartilage Marketa Kaucka 1,2 , Tomas Zikmund 3 , Marketa Tesarova 3 , Daniel Gyllborg 4 , Andreas Hellander 5 , Josef Jaros 6 , Jozef Kaiser 3 , Julian Petersen 2 , Bara Szarowska 2 , Phillip T Newton 1 , Vyacheslav Dyachuk 7 , Lei Li 1 , Hong Qian 8 , Anne-Sofie Johansson 8 , Yuji Mishina 9 , Joshua D Currie 10 , Elly M Tanaka 10 , Alek Erickson 11 , Andrew Dudley 11 , Hjalmar Brismar 12 , Paul Southam 13 , Enrico Coen 13 , Min Chen 14 , Lee S Weinstein 14 , Ales Hampl 6 , Ernest Arenas 4 , Andrei S Chagin 1,15 , Kaj Fried 7 *, Igor Adameyko 1,2 * 1 Department of Physiology and Pharmacology, Karolinska Institutet, Stockholm, Sweden; 2 Center for Brain Research, Medical University Vienna, Vienna, Austria; 3 Central European Institute of Technology, Brno University of Technology, Brno, Czech Republic; 4 Unit of Molecular Neurobiology, Department of Medical Biochemistry and Biophysics, Karolinska Institutet, Stockholm, Sweden; 5 Department of Information Technology, Uppsala University, Uppsala, Sweden; 6 Department of Histology and Embryology, Medical Faculty, Masaryk University, Brno, Czech Republic; 7 Department of Neuroscience, Karolinska Institutet, Stockholm, Sweden; 8 Department of Medicine, Karolinska Institutet, Stockholm, Sweden; 9 Department of Biologic and Materials Sciences, University of Michigan School of Dentistry, Ann Arbor, United States; 10 Center for Regenerative Therapies, Technische Universita ¨t Dresden, Dresden, Germany; 11 Department of Genetics, Cell Biology and Anatomy, University of Nebraska Medical Center, Omaha, United States; 12 Science for Life Laboratory, Royal Institute of Technology, Solna, Sweden; 13 John Innes Centre, Norwich, United Kingdom; 14 National Institute of Diabetes and Digestive and Kidney Diseases, National Institutes of Health, Bethesda, United States; 15 Institute for Regenerative Medicine, Sechenov First Moscow State Medical University, Moscow, Russia Abstract Cartilaginous structures are at the core of embryo growth and shaping before the bone forms. Here we report a novel principle of vertebrate cartilage growth that is based on introducing transversally-oriented clones into pre-existing cartilage. This mechanism of growth uncouples the lateral expansion of curved cartilaginous sheets from the control of cartilage thickness, a process which might be the evolutionary mechanism underlying adaptations of facial shape. In rod-shaped cartilage structures (Meckel, ribs and skeletal elements in developing limbs), the transverse integration of clonal columns determines the well-defined diameter and resulting rod-like morphology. We were able to alter cartilage shape by experimentally manipulating clonal geometries. Using in silico modeling, we discovered that anisotropic proliferation might explain cartilage bending and groove formation at the macro-scale. DOI: 10.7554/eLife.25902.001 Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 1 of 37 RESEARCH ARTICLE Introduction Cartilage is an essential skeletal and supportive tissue in our body. The shape and size of each cartilage element results from complex developmental processes; mesenchymal cells initially condensate, differentiate into chondrocytes, and then an orchestrated growth of the entire structure occurs (Goldring et al., 2006). Often, cartilage plays an important role as a developmental intermediate, such as during the endochondral growth of the long-bones (Mackie et al., 2008). Cartilage elements vary widely in their shapes: they may be simple shapes like rods or bars (Meckel, cartilage templates of the future long bones and ribs) or sheet-like structures (in the head), but can be extremely complicated with a huge number of irregular shapes (for instance, in the inner ear or pelvis). The geometrical properties of cartilage elements must be fine-tuned during the growth because cartilage provides indispensable structural support to the body during development. Yet, how this is achieved despite drastic changes in size is unclear. After early cartilage forms from mesenchymal condensations, growth typically occurs in all dimensions. However, the diversity of cell dynamics controlling precise early growth and shaping is not well studied. At the same time, the late growth of long rod-shaped cartilage elements in limbs is achieved through a mechanism of endochondral ossification that includes oriented cell dynamics in growth plate-like zones (Vortkamp et al., 1996). In the germinal zone of a growth plate, chondrocytes proliferate and produce progenies that form long streams oriented along the main axis of the forming skeletal element. Inside such streams, chondrocytes undergo flattening, oriented cell divisions and hypertrophy before dissipating and giving place to the forming bone (Nilsson et al., 2005), a process which is controlled by many signals (Kronenberg, 2003). This cell dynamic enables efficient extension of the skeletal element in a specific direction that coincides with the orientation of cell divisions in the proliferative zone (Abad et al., 2002). Growth plate disorders may result in dwarfism and other illnesses (De Luca, 2006). Some parts of the cartilaginous skull (e.g. the basisphenoid of the chondrocranium) also undergo endochondrial ossification in synchondroses, and significant growth of the cranial base is achieved through a similar mechanism (Hari et al., 2012;Wealthall and Herring, 2006). Synchondroses are mirror-image growth plates arising in the cranial base, which primarily facilitate growth in the anterio-posterior direction (Kettunen et al., 2006;Laurita et al., 2011; Nagayama et al., 2008;Young et al., 2006). Disorders in the development of synchondroses severely impact the elongation of the cranial base and often result in short-faced mutants and a general decrease of the cranial length (Ford-Hutchinson et al., 2007;Ma and Lozanoff, 1999). Insufficient or abnormal development of a cartilage element is one of the reasons for human craniofacial pathologies, providing a connection between the chondrocranium and facial bone geometry, size and placement (Wang et al., 1999). The growth mechanism operating in growth plates and synchondroses involves the transformation of the cartilage into the bone. Since growth plates or synchondroses are oriented towards a specific direction, the expansion of a cartilage in other dimensions is not clear from the mechanistic point of view and requires further investigation. For example, although it is well known that the mouse chondrocranium develops as 14 independent pairs of cartilage elements that form one united structure, the logic behind further shaping and scaling remains unclear (Hari et al., 2012). How these initially separated large cartilaginous elements form, grow and fine-tune their geometry, thickness and smoothness during development is still not completely understood. We hypothesized that accurate cartilage growth might require alternative cell dynamics that do not involve hypertrophy, ossification or growth plates. Such alternative cell dynamics may also contribute to the accuracy of scaling during cartilage growth. Scaling is a process of growth that maintains both the shape and the proportions of the overall structure. In nature, scaling often involves sophisticated principles of directional growth and a number of feedback mechanisms (Green et al., 2010). For instance, during bird development, the diversity in beak shape is constrained by the dynamics of proliferative zones in the anterior face (Fritz et al., 2014). Furthermore, scaling variations of beaks with the same basic shape result from signaling that controls the growth of the pre-nasal cartilage and the pre-maxillary bone (Mallarino et al., 2012). Indeed, in order to accurately scale a pre-shaped 3D-cartilaginous template both local isotropic and anisotropic cell dynamics may be required. Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 2 of 37 Research article Developmental Biology and Stem Cells To assess changes in the complete 3D anatomy of the face following cellular-level mechanistic studies we used a variety of approaches including micro-computed tomography (m-CT), genetic tracing with multicolor reporter mouse strains, multiple mutants and mathematical modelling. Most importantly, we reveal here how oriented clonal behavior in the chondrogenic lineage controls the overall geometry of the cartilage elements, and show that this geometry can be manipulated with molecular tools at various levels. Results Cartilage elements form and grow in all parts of the vertebrate body. The developing face provides a remarkable variety of cartilage geometries and sizes and, therefore, may serve as a sophisticated model system to study the induction of complex cartilaginous structures. The developing cartilaginous skull, the chondrocranium, displays a very complex geometry of mostly sheet-like cartilages that result from coordinated anisotropic growth in all dimensions. Such expansion of sheet-like cartilaginous tissue during embryonic development involves several mechanisms that were proposed in the past, including the formation and growth of cartilage at synchondroses, as well as at the apical growth zone. To understand the changes in dimensions of chondrocranium growth at major developmental stages, we took advantage of 3D reconstructions using m-CT enhanced with soft tissue contrasting (Figure 1). This approach enables the identification of various tissues and cell types in the embryo based on differential uptake of tungsten ions. We validated the m-CT visualization of embryonic cartilage by directly aligning stained histological sections with the 3D models (Figure 1—figure supplements 1 and 2). We analyzed the expansion of the chondrocranium due to synchondroses and found that despite a significant anterio-posterior elongation, synchondroses cannot entirely explain the growth dynamics in all directions: anterio-posterior, latero-medial and dorso-ventral vectors of growth (Figure 1A– C). Specifically, we found a complete absence of synchondroses and other endochondrial ossifications in the growing nasal capsule, even at the earliest postnatal stages, while membranous ossifications appeared well developed. The stereotypical clonal cell dynamics found in synchondroses (Figure 1D–I) did not appear during the development of the nasal capsule. Therefore, during the entire embryonic development, chondrocranium growth and shaping is largely aided by additional and unknown mechanisms of growth. To investigate another possible mechanism of growth, we examined the apical growth zone of the nasal capsule. To understand growth dynamics there, we birth-dated different regions of the facial cartilage using genetic tracing in Col2a1-CreERT2/R26Confetti and Sox10-CreERT2/R26Confetti embryos (Figure 2—figure supplement 1). Both Col2a1-CreERT2 and Sox10-CreERT2 lines recombine in committed chondrocyte progenitors and in mature chondrocytes. 3D analysis following tamoxifen injections at different developmental stages allowed us to identify the parts of the cartilage that develop from pre-existing chondrocytes and the regions generated from other cellular sources (Figure 2,Figure 2—figure supplement 1). As an example, after genetic recombination induced at E12.5, locations with high amount of traced cells show structures that come from preexisting cartilage, whereas areas comprising from non-traced cells present structures originating from de novo mesenchymal condensations. We discovered that important and relatively large geometrical features are produced from waves of fresh mesenchymal condensations induced directly adjacent to larger pre-laid cartilage elements between E13.5 and E17.5: this includes the frontal nasal cartilage, nasal concha, labyrinth of ethmoid and, consistent with previous suggestions, cribriform plate (Figure 2C and Figure 2—figure supplement 1A–K). These results cannot be safely inferred from 2D traditional histological atlases because of the complex geometry. Our results are complementary to the findings of McBratney-Owen and Morris-Key with coworkers, who demonstrated that the complete chondrocranium (including skull base) develops from 14 pairs of early independently induced large cartilaginous elements that fuse together during later development (McBratney-Owen et al., 2008). Here, we demonstrated how new adjacent mesenchymal condensation can increase the geometrical complexity of a single solid cartilaginous element. To substantiate our results, we took advantage of Ebf2-CreERT2/R26Tomato transgenic mouse line that can genetically label only a few selected patches of early mesenchyme in the cranial region. We wanted to test if some of these labelled mesenchymal patches can undergo chondrogenesis Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 3 of 37 Research article Developmental Biology and Stem Cells Figure 1. Visualizations of endochondrial ossifications in the chondrocranium during development. (A) 3D models of chondrocrania with visualized bone and hypertrophic cartilage. Note the absence of endochondrial ossifications in the nasal capsule between E14.5–17.5. Intramembraneous ossifications are not shown. (B) Width and length of the chondrocranium in E12.5–17.5 stages. (C) P2 stage model with visualized bone formation, hypertrophic zones and intramembraneous ossification in the nasal capsule. Clipping planes are applied for better visualizations of synchondroses. Figure 1 continued on next page Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 4 of 37 Research article Developmental Biology and Stem Cells independently and much later than most of the chondrocranium structure. If that would be the case, we could expect the formation of very sharp borders between the labelled and non-labelled cartilage due to the fusion of newly produced labelled cartilage with the old unlabeled one. If the local cartilage would form from labelled and unlabeled mesenchyme at the same time, the border would not form due to mesenchymal clone mixing that we observe when we label early neural crest. We injected Ebf2-CreERT2/R26Tomato animals with tamoxifen at E12.5 and analyzed the embryos at E17.5 (Figure 2—figure supplement 2). As a result, we discovered that the cartilage element connecting the inner ear with the basisphenoid was genetically traced, and demonstrated a very sharp border with non-traced cartilage (Figure 2—figure supplement 2C–D). m-CT data confirmed that this element develops entirely after E14.5 from newly formed mesenchymal condensations adjoining the chondrocranium (Figure 2—figure supplement 2A–B), and this might be related to differential regulation at the neural crest-mesodermal border (McBratney-Owen et al., 2008;Thompson et al., 2012). At the same time, the main structure of the chondrocranium is expanded in a very precise and symmetrical way due to unknown cellular and molecular mechanisms that cannot be explained by the freshly induced condensations, the apical growth zone, or even cell dynamics in synchondroses. Our m-CT results (Figure 2) show that various parts of the chondrocranium develop due to the growth of pre-existing cartilage not involving ossifications, while only additional features are induced in waves as de novo mesenchymal condensations that fuse with the main element during their maturation or expand in the process of ossification. We further focused on the developing nasal capsule because its growth does not involve synchondroses while the apical growth zone and adjoining mesenchymal condensation only partly provide for the growth and shaping modifications. The results obtained from comparisons of cartilaginous nasal capsules from different developmental stages showed that the shape of the structure is generally established by E14.5 (Figure 3, Figure 3—figure supplement 1,Figure 3—figure supplement 2,Video 1). Nevertheless, from E14.5 until E17.5 the cartilaginous nasal capsule is accurately scaled up with significant geometrical tuning (Figure 3A–B). Previous knowledge suggests that the underlying growth mechanism should be based on appositional growth of the cartilage during its transition to bone (Hayes et al., 2001; Li et al., 2017), however, numerous facial cartilages never ossify, but continue to grow. Tomographic reconstructions of sheet-shaped cartilage elements in the nasal capsule revealed extensive expansion of the cartilage surface area and overall volume (Figure 3E–F). Surprisingly, the thickness of the cartilaginous sheets did not change as much as the other dimensions during nasal capsule growth (Figure 3C–F,Figure 3—figure supplement 1,Video 1). Thus, the sheet-shaped cartilage expands mostly laterally (within the plane) during directional growth. Therefore, we expected that clonal analysis of the neural crest progeny (with Plp1-CreERT2/R26Confetti) and of early chondrocytes (with Col2a1-CreERT2/R26Confetti or Sox10-CreERT2/R26Confetti) would reveal clonal units (so called clonal envelopes) oriented longitudinally along the axis of the lateral expansion of the cartilage. Surprisingly, and contrary to this, clonal color-coding and genetic tracing demonstrated transversely oriented clones represented by mostly perpendicular cell columns or clusters formed by traced chondrocytes (Figure 4, especially A-C, Figure 4—figure supplement 1). To understand this process more in depth, we started with genetic tracing of the neural crest cells and their progeny in the facial cartilage with Plp1-CreERT2/R26Confetti (tamoxifen injected at E8.5). Clonal analysis and color-coding of neural crest-derived chondrogenic and non-chondrogenic ectomesenchyme showed intense mixing of neural crest-derived clones in any given location (Figure 4A– Figure 1 continued Corresponding raw CT data are presented in the lower part. (D–I) Clonal genetic tracing in synchondroses with Sox10CreERT2/R26Confetti; injected at E12.5 and collected at E17.5. (D) Schematic of synchondroses, (E) DAPI stained nuclei, (F–I) different clonal arrangements in various zones of progressing synchondroses. DOI: 10.7554/eLife.25902.002 The following figure supplements are available for figure 1: Figure supplement 1. Histological confirmation of m-CT results. DOI: 10.7554/eLife.25902.003 Figure supplement 2. Immunoand histological validation of cartilage contrasting obtained from m-CT analysis and subsequent 3D modelling. DOI: 10.7554/eLife.25902.004 Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 5 of 37 Research article Developmental Biology and Stem Cells Figure 2. Induction of initial shape and geometrical transformations of the facial chondrocranium. (A–C)m-CT-based 3D reconstruction of chondrogenic mesenchymal condensations and cartilage in the face of E12.5, E13.5, E14.5 and E17.5 embryos. (A) Mesenchymal condensations (yellow) segmented from E12.5 and E13.5 embryos and presented in frontal and top projections. Note that the basics of the facial chondrocranium are already established at the stage of mesenchymal condensations during the early development, while general geometry and fine details are tuned during further transformations. (B–C) Top, clipping plane + top and frontal projections of E14.5 (B) and E17.5 (C) developing facial chondrocranium. (C) Yellow color highlights the results of cartilage birth-dating experiments and shows the areas produced from de novo mesenchymal condensations that appear in successional waves after the primary cartilage (shown in green) is produced at previous stages (E14.5). Note that the shape of the facial chondrocranium develops as a result of incremental formation and additive fusion of new mesenchymal condensation with pre-existing cartilage. Red arrows indicate areas of cartilage which bend at later developmental stages (B,C) and red-outlined arrows indicate the same areas within the mesenchymal condensations at E13.5, prior to bending (A, bottom). DOI: 10.7554/eLife.25902.005 The following figure supplements are available for figure 2: Figure supplement 1. Genetic tracing serves as a tool for birth-dating of the cartilage during the embryonic development. DOI: 10.7554/eLife.25902.006 Figure 2 continued on next page Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 6 of 37 Research article Developmental Biology and Stem Cells C,Figure 4—figure supplement 1D–G) (Kaucka et al., 2016). At the same time, chondrogenic ectomesenchyme demonstrated the presence of transversely oriented doublets of genetically traced and also EdU-labeled cells already at E13.5 (Figure 4—figure supplement 1G (inserts) and H). Next, analysis of neural crest progeny in established cartilage highlighted the presence of perpendicularly oriented clonal doublets and columns (Figure 4A–C). Further analysis of EdU incorporation and genetic tracing with chondrocyte-specific Col2a1-CreERT2/R26Confetti and Sox10-CreERT2/ R26Confetti lines confirmed the existence of transversely oriented products of cell proliferation in the mature (E14.5-E17.5) cartilage (Figure 4D–F for EdU, Figure 4G,H–L and Figure 4—figure supplement 2 for lineage tracing). These results imply that cells in the sheet-shaped cartilage do not allocate daughter cells in lateral (longitudinal) dimensions as would be intuitively expected. Thus, simple lateral or unidirectional proliferation cannot account for the accurate scaling of the sheet-shaped cartilage in the face. Instead, the cartilage development from chondrogenic condensations is achieved by a cellular mechanism that involves intercalation of columnar clonal units. It was unclear to us why column-like structures, and no other shapes, are integrated into the sheet-shaped cartilage and how the fine surface is maintained during this mechanism of growth. To better understand possible mechanisms of accurate sheet-shaped cartilage surface development we modelled individual cell dynamics, in silico in 4D (3D + time) (Figure 5) (Hellander, 2015). We used this modelling to address two questions: firstly, under what conditions are clonal columns observed? Secondly, how is the sheet-like shape achieved by polarized or non-polarized cell divisions of singlecell thick layers and what are the controlling mechanisms? We tested a group of variables including: cell division speed, allocation of daughter cells in randomor defined directions, orientation cues in the tissue (equivalent to molecule gradients), as well as pushing/intercalating of the daughter cells during proliferation. We qualitatively compared the results from in silico simulations to our experimental clonal analysis from various genetic tracing experiments, in order to identify conditions in the model that were compatible with patterns observed in vivo. The results of the mathematical modelling suggested that the clonal dynamics observed in natural conditions requires polarity cues in the system, specifically, a two-sided gradient of signals would be required to precisely finetune cartilage thickness (Figure 5A–J). At the same time, some yet to be identified mechanism controls the average number of cell divisions in a column, further controlling columnar height and undoubtedly regulating the local thickness of the cartilage. Combined with the observed introduction of the transverse clonal columns, oriented cell proliferation can provide fine surface generation and scaling (Figure 4—figure supplement 2). Moreover, the model highlighted the elegance of cartilage design involving transverse columnar clones in the sheet-shaped elements: this logic enables the uncoupling of thickness control (depends on cell numbers within a clone) and lateral expansion (depends on the number of initiated clones), which are likely two molecularly unrelated processes in vivo. The absence of a gradient during in silico simulations led to the generation of 3D asymmetrical clusters instead of straight columns (even in conditions of highly synchronized cell divisions, and starting from a laterally space-constrained initial configuration - suggesting the promotion of vertical growth due to space-exclusion in the lateral direction) (Figure 5I,J). This, in turn, led to the formation of surface irregularities in the cartilage with subsequent loss of local flatness (heat-map diagram in Figure 5I,J). Importantly, lineage tracing also showed that for cartilaginous structures in the head with asymmetrical or complex irregular geometries, such as areas where several sheet-shaped cartilage elements were merged, clones were not constructed to perpendicular columns. In such locations, we identified irregular clonal clusters or randomly oriented clonal doublets, in accordance with the modelling results (Figure 6A–J). Thus, the shape and orientation of clones corresponds to the local geometry of the cartilage element. Next, we attempted to target a molecular mechanism that controls the flatness and sheet-like shape of the facial cartilaginous sheets. We discovered that activation of ACVR1 (BMP type one receptor, ALK2) in developing cartilage leads to a phenotype with targeted clonal micro-geometries Figure 2 continued Figure supplement 2. Genetic tracing based on Ebf2-CreERT2/R26Tomato serves as an indicator for structures developed from late mesenchymal condensations. DOI: 10.7554/eLife.25902.007 Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 7 of 37 Research article Developmental Biology and Stem Cells Figure 3. Facial chondrocranium undergoes major lateral expansion without extensive thickening during growth. (A) The 3D-model of E14.5 nasal capsule (blue) is placed onto the E17.5 model (green) for better presentation of growth-related changes. (B) Frontal clipping planes of 3D-models of nasal capsules at E14.5, E15.5, E16.5 and E17.5 (from left to right). Notice the mild changes in cartilage thickness as compared to the lateral expansion of the whole structure during growth. (C) Cartilage thickness heat-maps at E14.5 and E17.5 developmental stages. Less thick locations (color-coded in Figure 3 continued on next page Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 8 of 37 Research article Developmental Biology and Stem Cells (Figure 6K–P,R). We utilized a constitutively activated caALK2 transgene (Fukuda et al., 2006) together with genetic tracing in a way that every GFP-expressing cell is carrying constitutively active ACVR1. This experiment revealed a dramatic change of the shape of clonal envelopes, changing from straight perpendicular columns to disorganized spherical clusters inside the sheet-shaped cartilages of transgenic Sox10-CreERT2/R26caALK2-IRES-GFP embryos (Figure 6K–N). The ectopically activated ACVR1 resulted in the presence of clonal spherical clusters that interfered with the cartilage borders and caused the formation of ectopic bumps, swellings and other abnormal local shapes - in accordance with the mathematical modelling predictions (substantially resembling the condition with no gradient, see Figure 5I) (Figure 6O–P). All recombined cells in this caALK2 experiment became Sox9 + chondrocytes. There were no other cell types found to be GFP + , including perichondrial cells. This result indicates that BMP family ligands either produce the gradient that directs the orientated behavior of chondrocytes inside of the cartilage or, alternatively, that an experimental increase of BMP signaling renders the cells insensitive to the gradient established by other molecules. In any case, ACVR1 mutation can be used as a tool to change columnar arrangements into clusters (Figure 6N,R). The activation of ACVR1 by Sox10-CreERT2 starting from E12.5 occurred both in perichondrial cells and in chondrocytes (based on our genetic tracing results using Sox10CreERT2/R26Confetti). This later coincided with clonal bumps and bulging regions positioned mainly at the surface of sheet-shaped cartilaginous sheets (Figure 6L,N,P). These data also support the hypothesis that integration of clonal chondrocyte clusters into existing cartilaginous sheets likely depends on clonal shape and originates from the periphery of the cartilage. When this column-inserting process fails, the progeny of cells at the periphery of the cartilage forms ectopic bumps outside of the normal cartilage borders, and disrupts the flatness and straightness of cartilage surfaces. Next, we attempted to block the planar cell polarity (PCP) pathway to challenge the system and disrupt the formation of perpendicular columns in the flat or curved cartilaginous sheets. To do this we performed m-CT and EdU-incorporation analysis on Wnt/PCP mutants. Wnt/PCP pathway is well known for driving the cell and tissue polarity, and distinct facial phenotypes have appeared in Ror2, Vangl2 and Wnt5a homozygous mutants (Figure 7A). When EdU was administered 24 hr before embryo harvest, subsequent analysis showed no differences in the EdU-positive perpendicular clonal columns which formed within sheet-shaped facial cartilage of Figure 3 continued blue) correspond to intense growth zones shown in Figure 11. Dots show individual positions selected for precise measurements and demonstration on the graph shown in (D). Note that after E16.5 cartilage thickness remains relatively stable. (E–F) Cartilage surface area (E) and volume (F) expansion has been measured and compared between above mentioned stages. Note that there is a much greater increase in surface and volume (approximately 3-fold) than in thickness of the cartilage (less than 50%). DOI: 10.7554/eLife.25902.008 The following figure supplements are available for figure 3: Figure supplement 1. 3D models and wall thickness analysis of chondrocraniums at different developmental stages. DOI: 10.7554/eLife.25902.009 Figure supplement 2. Comparisons of the shape and size differences between developmental stages and Wnt/PCP mutants. DOI: 10.7554/eLife.25902.010 Video 1. 3D-models based on segmentation of mesenchymal condensations and mature cartilage from m-CT tomographic data. The first sequence illustrates wall thickness analysis results represented as a heatmap, starting from E12.5 (facial mesenchymal condensation) until E17.5 (facial cartilage). Cartilages and other soft tissues shrink during contrasting with phosphotungstic acid, and, thus, the reported metrics cannot be directly compared with biological samples treated in a different way. The following sequence shows facial chondrocranium models of Wnt/PCP mutants in comparison to the wild type. The last sequence shows the full chondrocranium at different embryonic stages, followed by 3D models of both the control embryo and Wnt5a mutant embryo at E17.5. DOI: 10.7554/eLife.25902.011 Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 9 of 37 Research article Developmental Biology and Stem Cells Figure 7. Deficiencies in Wnt/PCP pathway reshape the chondrocranium. (A)m-CT-based reconstructions of the facial chondrocranium of wild type control, Wnt5a -/- ,Vangl2 -/- ,Ror2 -/- and Vangl2 -/+ /Ror2 -/+ embryos at E15.5, with wall-thickness analysis (the row below). Clipping planes in the top projections show that all majorand fine structures (indicated by white arrows) are in place in the Wnt/PCP mutants. Red and blue arrows help to compare the width and the length of the chondrocranium. (B–D) Analysis of EdU incorporation in the facial sheet-shaped cartilage, 24 hr after the Figure 7 continued on next page Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 16 of 37 Research article Developmental Biology and Stem Cells Genetic tracing initiated during transition of condensations into cartilage resulted in clonal columns within both sheetand rod-shaped cartilage elements. This confirms that chondrogenic condensations undergo complex oriented cell dynamics during their development. Importantly, tracing of chondrocranium cartilage showed formation of transverse clonal columns as growth proceeded. Intercalation of newly born columns into pre-existing cartilage provided for the expansion potential in the sheet-shaped cartilage. This growth mechanism is very original and is not reported elsewhere so far. A few studies have demonstrated how clonal envelopes form in accordance with the general shape of the structure. These were mainly conducted on Drosophila imaginal wing disc or growing flower petals. In all cases the authors highlighted that the shape of clonal geometries correlates with the major vector of expansion in the growing structure (Green et al., 2010;Repiso et al., 2013; Strutt, 2005). This implies the presence of polarized activity that directs the shaping of the tissue. Here, we provided the first experimental evidence of how the control of the directional clonal expansion influences the shape of a vertebrate tissue on a large scale. Moreover, in the sheet-shaped cartilage the orientation of clonal domains, i.e. the columns, does not correspond to the vectors of major expansion, but rather serves for uncoupling lateral expansion control and thickness tuning. In line with that, the number of chondrocytes comprising the clonal column or cluster depends on Gsamediated signals. Variations in this number do not significantly affect the lateral dimensions of the whole sheet-shaped cartilage structure: the thickness of the cartilage becomes less while the general geometry and size stay preserved. Additionally, the shape and orientation of clonal envelopes in cartilage is partially controlled by BMP signaling, since micro-geometries of clones depend on activation of ACVR1. Based on these results, we assume that BMP ligands (because of cAlk2/ACVR1 phenotype affected clonal orientation) expressed around the regularly shaped cartilages may play a role similar to the in silico predicted gradients. Indeed, the expression of INHBA, BMP5 and BMP3 fit this expression profile quite well (according to Allen Developing Mouse Brain Atlas (http://developingmouse.brain-map.org) and Eurexpress (http://www.eurexpress.org) in situ public databases). At least, BMP5 is clearly expressed at the cartilage periphery and has been shown to affect the cartilage shape by David Kingsley lab (Guenther et al., 2008). Our experimental manipulations of planar cell polarity (PCP) pathway did not affect microgeometries and clonal domains, but strongly affected the chondrocranium shape on the macroscopical scale in several different ways. These phenotypes appeared to be rooted in pre-chondrogenic or early chondrogenic stages, and are based on distorted placement of mesenchymal condensations in the very early head. These experiments with Wnt/PCP mutants may potentially provide a better understanding of species-specific mechanisms of control and evolution of the facial shape on a macro scale. Regular shapes require regular cellular arrangements and clonal cell dynamics. It is not only sheet-shaped cartilage in the head that demonstrate geometric regularity; rod-shaped cartilage (Meckel, embryonic ribs and long cartilages in limbs) also has a regular shape. Regular clonal patterns, conceptually similar to those found in sheet-shaped cartilage, explain conservative tissue dynamics during formation and growth of cartilaginous rods. Indeed, genetic tracing experiments suggested that formation of clonal columns is important for the diameter control, while chondrogenic condensations at the very tip of the rod-shaped growing structures enable elongation. Similar to the cell dynamics in the sheet-shaped cartilage, this mechanism may provide for uncoupling of length versus diameter control. Such uncoupling may generally enable developmental and evolutional plasticity of cartilage size and shape. The mechanism controlling the thickness or diameter of sheet-shaped and rod-shaped cartilage elements not only includes spatially orientated behavior, but also involves the regulation of cell number within each chondrogenic clone. Immature chondrocytes are proliferatively active, while more Figure 7 continued pulse: control (B) and Wnt5a -/- mutant (C) embryos. Sparse clusters and columns of EdU + cells correspond to clonal arrangements previously shown with genetic tracing in the same locations. Note that the orientation of chondrocyte placement in the cartilage is not affected in the mutant embryos. Quantification is presented in a rose diagram in (D). For control, we refer to the rose diagrams in Figure 3. (E, F)m-CT-based 3D reconstruction of mesenchymal condensations at the developmental stage E12.5 in control (E) and Wnt5a -/- mutant (F) embryos shows their misplacement in a mutant. DOI: 10.7554/eLife.25902.017 Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 17 of 37 Research article Developmental Biology and Stem Cells Figure 8. Manipulation of cartilage thickness in GSamutant embryos. (A) Wall thickness was analysed in the m-CT segmented olfactory system of control (left) and GSa(G-protein subunit alpha) mutant (right). Large areas with decreased cartilage thickness are highlighted with a dashed line and white arrows. (B–G) Clonal genetic tracing of chondrocyte progenitors and chondrocytes induced by tamoxifen injection at E12.5 in Col2a1-CreERT2/ R26Confetti/GSa floxed/floxed embryos (C,E,G) and littermate controls (B,D,F) at E17.5. (H) Quantification of cartilage thickness in the olfactory system and Figure 8 continued on next page Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 18 of 37 Research article Developmental Biology and Stem Cells mature chondrocytes show decreased proliferation. Therefore, differentiation speed emerges as a concept which could regulate the organ shape by impinging on clone size, thereby altering the thickness or diameter of the cartilage. This concept is known to operate in the brain and other tissues with classical stem cell/transiently amplifying cell arrangements (Dı´az-Flores et al., 2006). Clonal genetic tracing and EdU labeling experiments suggested that the origin of clonal columns and clusters might be represented by the cells located at the periphery of forming cartilage. The spherical clusters of chondrocytes forming at the periphery of the cartilage in cALK2 mutant mice may suggest that the cell source is also located at the periphery and might be a perichondrial cell. Clonal relationships between perichondrial cells and columns of chondrocytes also support the hypothesis of perichondrial cells acting as a stem population during cartilage expansion. In general, the heterogeneity and multipotency of perichondrial cells is still unclear, although there are multiple studies showing the perichondrium as a source of chondrocytes and osteoblasts (Kobayashi et al., 2011;Li et al., 2017;Maes et al., 2010). In addition to this, the perichondrium might mediate non-autonomous effects in the cartilage in case of cAlk2 and GSaexperiments. Genetic tracing shows that some perichondrial cells always recombine with Sox10-, Plp1and Col2a1-CreERT2 lines, and, in case of functional experiments, may indirectly control some evens in more mature layers. Also, it is not clear how the fine border of the cartilage is set, and whether the perichondrial layer may play a key border-setting role during development and regeneration. This should be investigated further. Next, our results show that tuning of macro-geometries on a large scale can be achieved through a stage-specific placement of proliferative hot zones where new clonal domains intercalate into the main cartilage structure. Anisotropic heterogeneous proliferation is a powerful tool, which, together with polarity in the tissue and local patterning, can drive the organ shape development (Ben Amar and Jia, 2013;Campinho and Heisenberg, 2013). The localized growth zones provide for the general expansion and also bend the cartilage by creating local tensions that require mechanical relaxation and influence further development of the overall shape (Scho ¨tz et al., 2013). For probing such transformations of the sheet-shaped facial cartilage we applied an in silico model that was already successfully validated in a number of growth, shaping and scaling tasks (Green et al., 2010; Kennaway et al., 2011). Such a model was necessary to understand why the high and low proliferation zones are positioned in such a specific way. Indeed, the discovered distribution of proliferative zones in the whole nasal capsule did not help us per se with intuitive explanations of geometrical changes on the macro-scale. Despite this counter-intuitive dataset, the mathematical model provided an insight into the logic of the high and low proliferation zones in relation to a transition between investigated cartilage shapes. For example, it turned out that the position of lateral slow proliferation zones enables the generation of the symmetrical bends at the sides of the nasal capsule during transition from E13.5 to E14.5 developmental shapes. Furthermore, real material modelling confirmed the results predicted by the mathematical model, and generated lateral bends similarly to the native structure. The molecular mechanism controlling the dynamic distribution (patterning) of fast/slow proliferative zones in the cartilage is still unknown. It is likely linked to developmental signals from other tissues such as the olfactory epithelium or the mesenchyme surrounding the cartilage. Identification and validation of these signals will be essential in future studies and would involve a substantial combination of screening and functional approaches with transgenic animal models. Figure 8 continued basisphenoid from three independent experiments. Note the significant decrease of cartilage thickness in all analyzed locations. Oriented organization of the chondrocyte clones was not affected by GSaablation. The difference between control (mean = 5.9, sem = ±0.23, n = 4) and mutant (mean = 4.3, sem = ±0.25, n = 3) olfactory cartilage thickness is significant (p=0.0053). The difference between control (mean = 10.6, sem = ±0.83, n = 3) and mutant (mean = 5.7, sem = ±0.61, n = 3) basisphenoid cartilage thickness is significant (p=0.0087). Scale bars = 100 mm. (I) Graph showing that the regularity (straightness) of the cartilage correlates with the orientation of chondrocyte clones in the cartilage of Col2a1-CreERT2/R26Confetti/GSa floxed/floxed embryos. Angle acharacterizes the orientation of chondrocyte clones consisting of multiple cells as shown in a legend of a corresponding graph axis (y). Small values of acorrespond to highly oriented chondrocyte clones such as transverse clonal columns. Angle bis the angle between two opposite cartilage surfaces. Since sheet-shaped cartilage elements have almost parallel surfaces the angle bwas normally set between zero and 20 degrees. (J) GOM Inspect software was used to compare the shape of the nasal capsule between GSamutant and control embryo at E17.5. DOI: 10.7554/eLife.25902.018 Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 19 of 37 Research article Developmental Biology and Stem Cells The anisotropic proliferation can be an important evolutionary mechanism that is directly responsible for the differences in snout geometry in a variety of phylogenetic groups. Additionally, it might be important for understanding the development of the facial shape variation in humans (Sheehan and Nachman, 2014) as well as numerous pathologies (Afsharpaiman et al., 2013). One alternative way to fine-tune macro-geometry of a cartilage element is to continuously add on pre-shaped chondrogenic mesenchymal condensations from the pool of competent progenitors that Figure 9. Oriented cell dynamics during development of rod-shaped cartilage elements. (A) Genetic tracing in developing rib cartilage. Note the transverse pattern of chondrocyte clones. Dotted rectangle shows the area of magnified inset on bottom left. (B) 6 hr after EdU pulse in E14.5 embryo, transverse patterns were observed in ribs. Dotted rectangle shows the areas of magnified inset in bottom right corner. (C) Genetic tracing in developing rib cartilage shows transverse patterns. Dotted line represents area magnified in inset on the bottom right. (D–F) Genetic tracing (D,E) and 6 hr after EdU pulse (F) in developing digit cartilage of the upper limb. The areas of magnified insets (located on the right side) are outlined by dotted lines. (G) Transverse patterns in developing Meckel cartilage resulting from EdU pulse and analysis 6 hr after administration. (H) Genetic tracing shows transverse orientation of clonal chondrocytic columns in the Meckel cartilage. Dotted line shows the area magnified in the inset on the right. DOI: 10.7554/eLife.25902.019 Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 20 of 37 Research article Developmental Biology and Stem Cells Figure 10. Schematic overview of cartilage shaping and scaling processes. (1) Oriented cell divisions in the mesenchymal condensations give rise to the transverse columnar clones of chondrocytes. (2) Perichondrial cells may potentially give rise to chondrocytes. (3) Formation of new clonal columns and their integration into pre-existing cartilage leads to directed lateral expansion of the cartilage. The thickness of the sheet-shaped cartilage depends on the number of cells comprising the column, while the lateral expansion depends on the number of clonal columnar units engaged. (4) Geometry of the Figure 10 continued on next page Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 21 of 37 Research article Developmental Biology and Stem Cells are retained until late developmental stages. As we demonstrated, the formation of adjoining mesenchymal condensations occurs in sheet-shaped cranial cartilage. In the developing face, new chondrogenic condensations are responsible for introducing geometrically complicated fine details. Such mechanisms may also operate during amphibian metamorphosis, when most of the postmetamorphic cranial cartilage develops de novo and not from the pre-metamorphic cartilaginous elements (Kerney et al., 2012). Taken together, we discovered important novel principles explaining the growth and shaping of cartilaginous structures. Further studies should focus, amongst other things, on the soluble signals emanating from other embryonic structures which influence the oriented behavior or proliferation of chondrogenic clones. Materials and methods Mouse strains and animal information All animal (mouse) work has been approved and permitted by the Ethical Committee on Animal Experiments (Norra Djurfo ¨rso ¨ksetiska Na ¨md, ethical permit N226/15 and N5/14) and conducted according to The Swedish Animal Agency´s Provisions and Guidelines for Animal Experimentation recommendations. Genetic tracing mouse strains Plp1-CreERT2 (RRID:MGI:4837112) and Sox10CreERT2 were previously described (Laranjeira et al., 2011;Leone et al., 2003;Yu et al., 2013). Plp1-creERT2,Sox10-creERT2 and Col2a1-CreERT2 (RRID:IMSR_JAX:006774) (Nakamura et al., 2006) (obtained from laboratory of S. Mackem, NIH) strains were coupled to R26Confetti (RRID: IMSR_JAX:017492) mice that were received from the laboratory of Professor H. Clevers (Snippert et al., 2010). The Stop flowed/floxed caAlk2-IRES-GFP strain from the laboratory of Y. Mishina (Fukuda et al., 2006) was coupled to Sox10-CreERT2. The Ebf2-CreERT2 (RRID:MGI:4421811) strain was obtained from the laboratory of H. Qian, KI, and was coupled to R26Tomato. The Gsa floxed/floxed strain was obtained from the laboratory of L. Weinstein (Sakamoto et al., 2005). Female mice which were homozygous for the reporter allele [Gt(ROSA)26Sortm4(ACTB-tdTomato,-EGFP)Luo/J; Jackson Laboratories] (Muzumdar et al., 2007) were coupled to homozygous Col2a1::creER T males [FVB-Tg (Col2a1-cre/ERT)KA3Smac/J; Jackson Laboratories] (Feil et al., 1997;Nakamura et al., 2006). To induce genetic recombination to adequate efficiency, pregnant females were injected intraperitoneally with tamoxifen (Sigma Aldrich, St.Louis, MO, T5648) dissolved in corn oil (Sigma Aldrich, C8267). Tamoxifen concentration ranged from 1.5 to 5.0 mg per animal in order to obtain a range of recombination efficiencies. Wnt5a, Vangl2 and Ror2 full knock-out embryos were obtained from heterozygous parents (Gao et al., 2011;Yamaguchi et al., 1999) at the expected Mendelian proportions. Immunohistochemistry For embryo analyses, heterozygous mice of the relevant genotype were mated overnight, and noon of the day of plug detection was considered E0.5. Mice were sacrificed with isoflurane (Baxter, Deerfield, IL, KDG9623) overdose, and embryos were dissected out and collected into ice-cold PBS. Subsequently, the samples were placed into freshly prepared 4% paraformaldehyde (PFA) and depending on the developmental stage they were fixed for 3–6 hr at +4˚C on a roller. Embryos were subsequently cryopreserved in 30% sucrose (VWR, Radnor, PA, C27480) overnight at +4˚C, embedded in OCT media (HistoLab, Serbia, 45830) and sections cut of between 14 mm to 200 mm on a cryostat (Microm International, Germany), depending on the following application. If needed, sections were stored at 20˚C after drying for 1 hr at room temperature, or processed immediately Figure 10 continued clonal unit corresponds to the overall macro-geometry of the cartilage. Regular clonal units correspond to regular shapes of the cartilage. (5) Chondrogenic mesenchymal condensations are sequentially induced to provide fine details and shape modifications during chondrocranium growth. Upon their maturation, they fuse with the main structure of the chondrocranium. Anisotropic proliferation and specifically positioned proliferative zones further assist the shaping process by imposing physical tensions and curves. (6) Rod-shaped cartilage elements also show the regular clonal patterns that result from the transverse orientation of cell divisions and daughter cell allocations that account for the diameter regulation. DOI: 10.7554/eLife.25902.020 Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 22 of 37 Research article Developmental Biology and Stem Cells after sectioning. Primary antibodies used were: goat anti-GFP (FITC) (Abcam, UK, 1:500, RRID:AB_ 305635), rabbit anti-Sox9 (Sigma Aldrich, 1:1000, RRID:AB_1080067), rabbit anti-Sox5 (Abcam, 1:500, RRID:AB_10859923), sheep anti-ErbB3 (RnD Systems, Minneapolis, MN, 1:500, RRID:AB_ 2099728). For detection of the above-mentioned primary antibodies we utilized 405, 488, 555 or 647-conjugated Alexa-fluor secondary antibodies produced in donkey (Invitrogen, Carlsbad, CA, 1:1000, RRID:AB_162543, RRID:AB_141788, RRID:AB_141708, RRID:AB_142672, RRID:AB_2536183, RRID:AB_141844,). Sections were mounted with 87% glycerol mounting media Figure 11. Analysis of proliferation identifies specific proliferative regions in nasal capsule. (A–B) Analysis of EdU incorporation 24 hr after the pulse on a transversal section of the facial chondrocranium at E14.5. Notice the distinct proliferative zones in the cartilage that correlate with intense EdU labelling in perichondrial locations shown by arrowheads in (B). (C) Mapping of distinct growth zones onto 3D models of mesenchymal condensations (E13.5) and cartilage (E14.5-E15.5) in the developing face. (D–F) Frontal transversal sections at different developmental stages include proliferative zones within the chondrocranium with EdU incorporation. Scale bars = 100 mm. DOI: 10.7554/eLife.25902.021 Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 23 of 37 Research article Developmental Biology and Stem Cells Figure 12. Modeling of uneven growth in the shaping of the nasal capsule. (A) In silico geometrical transformations of the nasal capsule-like anlage at E13.5 following various scenarios including: anisotropic oriented growth (following polarization introduced by the anterio-posterior gradient shown in green), non-polarized isotropic growth (no anterio-posterior gradient), the presence of fixed midline (simulation of septum and central groove), condition with the unfixed midline (only central groove), conditions with or without slowly growing lateral regions (shown in purple). Note that in Figure 12 continued on next page Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 24 of 37 Research article Developmental Biology and Stem Cells (Merck, Germany) or in Vectashield Antifade Mounting Medium with DAPI (Vector Laboratories, Burlingame, CA, RRID:AB_2336790). EdU incorporation analysis EdU (Life Technologies, Carlsbad, CA) was injected intraperitoneally into the pregnant females (65 mg per gram of body mass) either 6or 24 hr before the embryos were harvested. Cells with incorporated EdU were visualized using Click-iT EdU Alexa Fluor 647 Imaging Kit (Life Technologies) according to the manufacturer’s instructions. Microscopy, volume rendering, image analysis and quantifications Confocal microscopy was performed using Zeiss LSM710 CLSM, Zeiss LSM780 CLSM and Zeiss LSM880Airyscan CLSM instruments. The settings for the imaging of Confetti fluorescent proteins were previously described (Snippert et al., 2010). The imaging of the confocal stack was done with a Zeiss LSM780 CLSM, Plan-Apochromat 3 10x/0.45 M27 Zeiss air objective. Histological staining Slides were stained for mineral deposition using von Kossa calcium staining: 5% silver nitrate solution was added to the sections at a room temperature and exposed to strong light for 30 min. After that the silver nitrate solution was removed, and slides were washed with distilled water for three times during 2 min. 2.5% sodium thiosulphate solution (w/v) was added to the sections and incubated for 5 min. Slides were again rinsed for three times during 2 min in distilled water. The sections were then counterstained using Alcian blue. Alcian blue solution (0.1% alcian blue 8GX (w/v) in 0.1 M HCl) was added to the tissue for 3 min at room temperature and then rinsed for three times during 2 min in distilled water. Slides were then transferred rapidly into incrementally increasing ethanol concentrations (20%, 40%, 80%, 100%) and incubated in 100% ethanol for 2 min. Finally, the slides were incubated in two xylene baths (for 2 min and then for 5 min) before mounting and analysis. Figure 12 continued condition with polarized anterio-posterior growth the anterior elongation of the structure is more prominent and faster as compared to the condition with non-polarized isotropic growth. The lateral bends are induced by slow proliferating lateral regions. In the center and on the right, the real nasal capsules are shown with mapped fast and slow growing regions. (B) In the condition with isotropic growth and introduced slowly growing regions, we observe the formation of lateral bends (red line) analogous to the lateral bends in the real nasal capsule at E14.5 (shown on the left). In this condition the midline is fixed, and the ventral groove forms straight. (C) In conditions with no fixed midline we observe the formation of the central groove, correct bending of the central groove (red line) and overall flattening of the simulated structure similar to the real object (on the left). (D) Simulation with no midline and central groove. Note the inverted bend (red line) and the absence of the correct flattening of the structure. Despite the absence of the midline, the lateral bends are successfully induced by the slow growing regions (purple), analogous to the real nasal capsule. (E) Material elastic modelling shows how the third dimension (bending) emerges from changes and tensions in plain 2D structure during imitated anisotropic growth. (E, left panel) Initial modelling conditions: completely flat X-ray film with the cut slot in the middle for fitting the imitated flat growth zone, which is also made from X-ray film. (E, right panel) When the growth zone is inserted into the slit, the whole structure bends to accommodate the tensions. (F) Real material (plastic film)-based simulation of isotropic growth was based on uneven shrinking during intense heating. Black painted regions uptake heat more efficiently and shrink faster. The attached edges of the shrinking zone cause bending of the entire structure. Two lateral black stripes were painted on top of the trapezoid as an analog to lateral slowly proliferating zones in nasal capsule. Note the similarity of resulting bends to the lateral bends in real nasal capsule at E14.5. DOI: 10.7554/eLife.25902.022 Video 2. Simulations of shape transitions of the nasal capsule-like 3D object under different growth conditions. Notice the formation of the lateral bends corresponding to the real nasal capsule shape development from E13.5 to E14.5 occur only in the condition with slow growing purple zones. These bends form independently from isotropic or anisotropic modality of growth. DOI: 10.7554/eLife.25902.023 Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 25 of 37 Research article Developmental Biology and Stem Cells Nagayama M, Iwamoto M, Hargett A, Kamiya N, Tamamura Y, Young B, Morrison T, Takeuchi H, Pacifici M, Enomoto-Iwamoto M, Koyama E. 2008. Wnt/beta-catenin signaling regulates cranial base development and growth. Journal of Dental Research 87:244–249. doi: 10.1177/154405910808700309,PMID: 18296608 Nakamura E, Nguyen MT, Mackem S. 2006. Kinetics of tamoxifen-regulated cre activity in mice using a cartilagespecific CreER(T) to assay temporal activity windows along the proximodistal limb skeleton. Developmental Dynamics 235:2603–2612. doi: 10.1002/dvdy.20892,PMID: 16894608 Nilsson O, Marino R, De Luca F, Phillip M, Baron J. 2005. Endocrine regulation of the growth plate. Hormone Research in Paediatrics 64:157–165. doi: 10.1159/000088791,PMID: 16205094 Radszuweit M, Block M, Hengstler JG, Scho ¨ll E, Drasdo D. 2009. Comparing the growth kinetics of cell populations in two and three dimensions. Physical Review E 79:051907. doi: 10.1103/PhysRevE.79.051907 Repiso A, Bergantin ˜os C, Serras F. 2013. Cell fate respecification and cell division orientation drive intercalary regeneration in Drosophila wing discs. Development 140:3541–3551. doi: 10.1242/dev.095760,PMID: 239031 86 Sakamoto A, Chen M, Kobayashi T, Kronenberg HM, Weinstein LS. 2005. Chondrocyte-specific knockout of the G protein G(s)alpha leads to epiphyseal and growth plate abnormalities and ectopic chondrocyte formation. Journal of Bone and Mineral Research 20:663–671. doi: 10.1359/JBMR.041210,PMID: 15765186 Scho ¨tz EM, Lanio M, Talbot JA, Manning ML. 2013. Glassy dynamics in three-dimensional embryonic tissues. Journal of the Royal Society Interface 10:20130726. doi: 10.1098/rsif.2013.0726,PMID: 24068179 Sheehan MJ, Nachman MW. 2014. Morphological and population genomic evidence that human faces have evolved to signal individual identity. Nature Communications 5:4800. doi: 10.1038/ncomms5800,PMID: 252262 82 Snippert HJ, van der Flier LG, Sato T, van Es JH, van den Born M, Kroon-Veenboer C, Barker N, Klein AM, van Rheenen J, Simons BD, Clevers H. 2010. Intestinal crypt homeostasis results from neutral competition between symmetrically dividing Lgr5 stem cells. Cell 143:134–144. doi: 10.1016/j.cell.2010.09.016,PMID: 20887898 Strutt D. 2005. Organ shape: controlling oriented cell division. Current Biology 15:R758–R759. doi: 10.1016/j. cub.2005.08.053,PMID: 16169474 Tesar ˇova ´M, Zikmund T, Kaucka´ M, Adameyko I, Jarosˇ J, Palousˇek D, S ˇkaroupka D, Kaiser J. 2016. Use of micro computed-tomography and 3D printing for reverse engineering of mouse embryo nasal capsule. Journal of Instrumentation 11:C03006. doi: 10.1088/1748-0221/11/03/C03006 Thompson H, Ohazama A, Sharpe PT, Tucker AS. 2012. The origin of the stapes and relationship to the otic capsule and oval window. Developmental Dynamics 241:1396–1404. doi: 10.1002/dvdy.23831,PMID: 2277 8034 Van Liedekerke P, Palm MM, Jagiella N, Drasdo D. 2015. Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results. Computational Particle Mechanics 2:401–444. doi: 10.1007/ s40571-015-0082-3 Vortkamp A, Lee K, Lanske B, Segre GV, Kronenberg HM, Tabin CJ. 1996. Regulation of rate of cartilage differentiation by indian hedgehog and PTH-related protein. Science 273:613–622. doi: 10.1126/science.273. 5275.613,PMID: 8662546 Wang Y, Spatz MK, Kannan K, Hayk H, Avivi A, Gorivodsky M, Pines M, Yayon A, Lonai P, Givol D. 1999. A mouse model for achondroplasia produced by targeting fibroblast growth factor receptor 3. PNAS 96:4455– 4460. doi: 10.1073/pnas.96.8.4455,PMID: 10200283 Wealthall RJ, Herring SW. 2006. Endochondral ossification of the mouse nasal septum. The Anatomical Record Part A: Discoveries in Molecular, Cellular, and Evolutionary Biology 288:1163–1172. doi: 10.1002/ar.a.20385, PMID: 17031811 Yamaguchi TP, Bradley A, McMahon AP, Jones S. 1999. A Wnt5a pathway underlies outgrowth of multiple structures in the vertebrate embryo. Development 126:1211–1223. PMID: 10021340 Young B, Minugh-Purvis N, Shimo T, St-Jacques B, Iwamoto M, Enomoto-Iwamoto M, Koyama E, Pacifici M. 2006. Indian and sonic hedgehogs regulate synchondrosis growth plate and cranial base development and function. Developmental Biology 299:272–282. doi: 10.1016/j.ydbio.2006.07.028,PMID: 16935278 Yu K, McGlynn S, Matise MP. 2013. Floor plate-derived sonic hedgehog regulates glial and ependymal cell fates in the developing spinal cord. Development 140:1594–1604. doi: 10.1242/dev.090845,PMID: 23482494 Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 32 of 37 Research article Developmental Biology and Stem Cells Appendix 1 Individual based model (IBM) for cartilage dynamics In order to model and illustrate the growth of the cartilage on the cellular level, we developed an individual-based model incorporating cell proliferation (including displacement of surrounding cells via pushing) (Hellander, 2015). The model is stochastic because we want to be able to capture effects of e.g. synchronicity in cell division and the degree of determinism needed to achieve and ordered columnar growth of the structure. We do this by letting the cell division times, direction of allocation of daughter cells after cell division, etc., be random variables. This document describes the details of the model and its implementation. There are several popular modeling frameworks for simulating interacting cells. In the cellular potts model (Graner and Glazier, 1992) a single biological cell can be composed of multiple lattice sites making it possible to use a more detailed description of cell shape and include more detailed description of more mechanical properties. In off-lattice center based models cells are often modeled as spheres with pair-wise interactions and a force-based description to evolve the system dynamics, for an overview see (Van Liedekerke et al., 2015). Vertex models can offer even more realistic models of cell mechanics (Fletcher et al., 2013) but become expensive and complicated in three space dimensions. Rather than these more comprehensive mechanical models, we use a simplistic rule-based, on-lattice stochastic cellular automaton (CA). In the language of a recent review (Van Liedekerke et al., 2015) our model falls into the category of a Type B CA. These types of models are widely used in e.g. cancer tumor modeling and for simulation of monolayers and spheroids (Radszuweit et al., 2009). The simulation code is written in Python, relies on the PyURDME package for spatial stochastic simulations (www.pyurdme.org) and it is freely available for download from www.github.com/ahellander/multicell (Hellander, 2015) under the GPLv3 license. A copy is archived at https://github.com/elifesciences-publications/ multicell. The basic entities in our simulation are Agents and Events. An Agent is a model (implemented as a Python class) of an individual (cell). Events simulate discrete state changes involving one or several agents. They occur at a certain time (assuming no other event involving the same agents occurs first). They rely on rules that specify how and under what conditions the event is to be executed. A simulation is initialized by creating the initial population of agents and events, and then creating a priority queue (in our case implemented using a heap data structure using the Python module’ heapq’). In each iteration of the algorithm, the event with the shortest time is popped from the queue and executed (assuming that all of its rules and conditions can be satisfied), the system time updated, new events derived from any newly created agents are created and inserted into the queue, and all existing events affected by changes in the agents or the system state are updated. Each agent occupies one voxel of a tessellation of 2D or 3D space, and each lattice site can only accommodate one agent. Following the recommendation in (Radszuweit et al., 2009) a Delaunay triangulation is used. The mesh resolution is chosen such that the average voxel size is close to the desired cell size (~7mm radius) taken from the experiments. Being a lattice model, the shape and volume of the cells are a lattice property and are given by the dual grid (Voronoi cells in the case of a Delaunay triangulation). This is illustrated in Appendix 1—figure 1A. The interpretation is that individual agents occupy the dual elements (dashed lines) of an unstructured triangular (2D) or tetrahedral primal mesh (3D) (solid lines). Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 33 of 37 Research article Developmental Biology and Stem Cells Appendix 1—figure 1. (A) Individual cells are modeled by a number of properties such as their color and distributions for cell division times. The positions of cells in space are tracked on an underlying unstructured lattice, or grid. The edges in the primary mesh (solid lines) connect vertices (black dots). A biological cell is modeled by the volume made up of the dual elements (dashed lines), connecting triangle (2D) or tetrahedral (3D) centers and edge or face centers. For visualization purposes, in 3D space, we plot cell individuals as colored spheres with radius equal to the sphere with equal volume as the dual element. (B) Size distribution for the mesh elements for the geometry and mesh used in the simulations in the Figure 5. DOI: 10.7554/eLife.25902.025 The individual agents – colored coded cells Each individual cell is modeled as an individual agent with the following properties: .Color (a label used to track the lineage). .Mean cell division time, p. .Variance in cell division time, s2 p: When visualized in 2D, we draw cells as polygons (the actual dual cells) and in 3D for practical reasons we visualize them as spheres centered on the vertices of the primal mesh, with radius chosen such that the volume corresponds to the volume of the corresponding dual element. On the unstructured mesh, there will be a size distribution for the mesh elements, i.e. there will be a small variation in the size associated with each lattice site, see Appendix 1— figure 1B. Cell proliferation Cell division time The time until a cell, or individual, divides, is assumed to be a random variable. Although a multi-stage model of cell division can give rise to an an Erlang distribution (Radszuweit et al., 2009) which is found to match experimental data for another system, we are not calibrating our model to experimental data on cell division time distributions. The dividing cell (referred to as the mother cell) create a daughter cell after a normally distributed waiting time tD. tD~Nd;s2 d   (1) Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 34 of 37 Research article Developmental Biology and Stem Cells where the mean division time mand the variance s2are parameters of the model to be supplied as input to the simulation. As a measure of the degree of variability in cell division times, we use the standard deviation over the mean, f¼sd d (2) The smaller value of f, the more deterministic and synchronized the cell cycles of individual cells are. The way our model is set up, there are no events that leads to the recalculation of a cell’s division time. A cell gets assigned a time to division at creation and each cell then divide according to its internal clock irrespective of if it gets pushed etc. After division, the daughter cell needs to be deposited on the grid. The division direction, or receiving voxel, is sampled according to a discrete distribution. In the simplest case, all directions of division are equally probable and the direction distribution is uniform. In the general case, weights are assigned according to an external, deterministic gradient. Cell division direction The division direction is also a random variable and the number of possible directions are given by the connections to the neighbors on the grid. Each individual has a property that sets its polarization, represented by a normalized vector ppointing in the preferred direction of division. With no polarization, each possible division direction is equally probable. With polarization, the probability to divide in a certain direction is biased by the gradient. The weights for sampling the division direction are taken to be wij ¼dij maxjdij   !b ;(3) dij ¼g xj  g xi ð Þ hij (4) where xiis the position of the vertex in the grid for which the agent resides, and hij is the length of the edge connecting grid points xiand xjand g xð Þ is a given concentration profile. The parameter b0dictates how perfectly the cells become polarized by the concentration profile g xð Þ. A value b¼0leads to equal probabilities for all directions, and very large value of bmeans that the division direction will always be in the direction of the maximal value of the gradient (the division direction becomes deterministic in the direction of the maximal gradient in the limit b!¥. Values in between the extremes describes an increasing precision in polarization axis alignment with the gradient field. Cell pushing If the receiving lattice for the daughter cell site is empty (i.e. occupied by matrix), it is simply deposited there. When the daughter cell cannot be placed on a free lattice site, there is an attempt to reorganize the structure by pushing neighboring cells to make room for it. The procedure is illustrated in Appendix 1—figure 2. The probability for the displaced cell to move to a given neighboring grid point depends on the direction of pushing. Let emd be the vector along the edge connecting the mother cell Cmand the daughter cell Cd, pointing towards the daughter cell. Let edk be the unit vector along the edge connecting the daughter cell Cdand one of its neighbors, Ck. The weight for moving the displaced cell to the neighbor with index nis given by Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 35 of 37 Research article Developmental Biology and Stem Cells Appendix 1—figure 2. When a cell divides (yellow cell), the daughter cell will push surrounding cells (blue) to make room for the progeny (A). The direction of pushing (and what cell to push) is determined by a combination of the directivity of the original division or pushing direction and of a penalty for pushing an occupied site. The penalty is governed by a parameter c that dictates how much to favor pushing into a free lattice site. In the figure, pushing the blue cell at location Cdto Clwill be favored over Ckbased on the penalty if Ck>0, but the site Ck if favored based on the directionality of the push. If Ckis selected by the probabilistic algorithm, it will in turn complete a pushing event, leading to a pushing chain that continues until a cell is pushed into an empty site. After such a pushing chain has been completed, the site sampled for the daughter cell will be free, and the newly created daughter cell can be inserted (B). DOI: 10.7554/eLife.25902.026 wk¼max sk1cIn ð Þ;0ð Þ (5) where sk¼emd edk (6) That is, to account for the directionality implied by the originally sampled division direction emd, the probability of the pushed cell to move to an adjacent lattice site ckis proportional to the scalar projection of the vector connecting the pushed cell and the vertex at index k onto the direction vector of the push. cis a penalty parameter in the interval [0, 1] that dictates the degree of resistance in pushing the cell into an already occupied site. If c¼1, it is not possible to push a cell into an occupied site, and if c¼0, there is no resistance what so ever and only the direction of the push affects the displacement direction. Any value in between is a tradeoff between the two extremes. In Appendix 1—figure 2, for example, there is a high probability to push Cjto Ck, due to edk being almost parallel to emd (high sk), but depending on the value of c, the site Ckmay be sampled instead since that site is empty, even if the directionality contributions is smaller ðsl<skÞ. Once a neighbor has been displaced, the pushed cell moves into that lattice site and the daughter cell gets deposited on the now free lattice site of the displaced cell. The displaced cell, may in turn then go on to displace additional cells and this procedure is repeated until a cell gets pushed into an unoccupied site. Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 36 of 37 Research article Developmental Biology and Stem Cells Measure of order in the cartilage model We are interested in assessing what factors are the main determinants to the degree of ordered columnar growth in the cartilage sheet patches. To that end, we postulate that a perfectly ordered structure consists of clonal columns growing straight and directed along the axis perpendicular to the initial condition starting plane (Figure 5G). We then use the following metric to quantify the degree of order in the structure S¼1 CX C j¼1 1 NjX Nj i¼1 sij (7) with sij sij ¼ jvi.yj(8) where yis a unit vector perpendicular to the initial condition plane, and viare normalized vectors joining two consecutive points (sorted by y-coordinate) in clones with the same color. Njis the number of cells of a given color minus one (the number of vectors), and Cis the number of unique clones tracked. This is illustrated graphically in Figure 5A. With this metric, a score of S¼1would mean that all columns are perfectly aligned to the main growth axis and a score of S¼0would mean that they are all perpendicular to it. We use this metric to score realizations of the process either in the absence of a gradient, or when the gradient is uniform in planes parallel to the center plane, so that in the case of a perfect polarization, cells should all deposit their daughter cells perpendicular to the center plane. Kaucka et al. eLife 2017;6:e25902. DOI: 10.7554/eLife.25902 37 of 37 Research article Developmental Biology and Stem Cells