Shaping of a three-dimensional carnivorous trap through modulation of a planar growth mechanism

Leaves display a remarkable range of forms, from flat sheets with simple outlines to cup-shaped traps. Although much progress has been made in understanding the mechanisms of planar leaf development, it is unclear whether similar or distinctive mechanisms underlie shape transformations during development of more complex curved forms. Here, we use 3D imaging and cellular and clonal analysis, combined with computational modelling, to analyse the development of cup-shaped traps of the carnivorous plant Utricularia gibba. We show that the transformation from a near-spherical form at early developmental stages to an oblate spheroid with a straightened ventral midline in the mature form can be accounted for by spatial variations in rates and orientations of growth. Different hypotheses regarding spatiotemporal control predict distinct patterns of cell shape and size, which were tested experimentally by quantifying cellular and clonal anisotropy. We propose that orientations of growth are specified by a proximodistal polarity field, similar to that hypothesised to account for Arabidopsis leaf development, except that in Utricularia, the field propagates through a highly curved tissue sheet. Independent evidence for the polarity field is provided by the orientation of glandular hairs on the inner surface of the trap. Taken together, our results show that morphogenesis of complex 3D leaf shapes can be accounted for by similar mechanisms to those for planar leaves, suggesting that simple modulations of a common growth framework underlie the shaping of a diverse range of morphologies.


Introduction
Many plant and animal organs, such as leaves, flowers, hearts, and wings, derive from tissue sheets. A general question in developmental and evolutionary biology is how tissue sheets are shaped to create such a diversity of forms. A good illustration is leaf development. Leaves exhibit remarkable variation in shape, from simple or compound planar forms to convoluted three-dimensional forms such as those of pitcher plants. The molecular genetic control of leaf shape has been extensively studied for planar forms, with key genes modifying leaf shape identified [1][2][3][4][5][6]. Variation in patterns of gene activity has also been shown to underlie variation in leaf shape between species [7][8][9][10]. Clonal analysis and tracking and monitoring cell division have further revealed spatiotemporal variation in patterns of division and growth, and led to the formulation of models for how shape arises through local variations in rates and orientations of growth [11,12]. However, it is unclear how these models for planar leaf development are related to morphogenetic changes in highly curved 3D leaf forms, such as epiascidiate (cup or tubular-shaped) leaves.
Epiascidiate leaves have evolved four times independently: in the families Nepenthaceae, Sarraceniaceae, Cephalotaceae, and Lentibulariaceae [13,14]. In all these cases, the epiascidiate form is associated with nutrient acquisition from animals (carnivory). Based on comparative anatomy, the inner surface of the epiascidiate leaf is believed to be equivalent to the adaxial surface of a planar leaf, whereas the outer surface is equivalent to the abaxial surface [4,13,15,16]. The petiole of the epiascidiate leaf inserts on the abaxial side, similar to the situation for a peltate leaf [13]. However, the mechanism by which the epiascidiate leaf is initially formed and then shaped during development is poorly understood. Here, we address the developmental mechanisms controlling the second aspect, involving shaping of a highly curved sheet.
In the genus Utricularia (Lentibulariaceae), epiascidiate leaves, termed traps, use suction to catch prey, requiring highly coordinated morphogenesis to ensure that the opening and closing mechanisms operate effectively. Utricularia has several advantages for analysis of epiascidiate leaf development [17]. The traps are transparent and only a few millimetres long, making them convenient for imaging. Much of the trap comprises only two cell layers [18,19], compared to approximately seven cells for Arabidopsis leaves, simplifying growth analysis. The genome of Utricularia gibba is among the smallest in plants (100 Mb) and has been fully sequenced, providing a resource for molecular genetic and evolutionary studies [20][21][22][23][24][25][26][27]. Utricularia is also a large genus, comprising about 235 species with varying trap shapes, allowing for comparative analysis [28][29][30][31].
Snapshots, scanning electron micrographs, and drawings of Utricularia traps at various developmental phases have been described [18,[32][33][34][35]. However, quantitative growth and cellular analysis of morphogenesis have not been carried out. Such studies require the development of transformation methods for introducing fluorescent proteins to mark cell membranes or clones in Utricularia, followed by 3D imaging at different developmental stages. Moreover, models need to be developed for how changes in three-dimensional shape and curvature arise, and predictions of these models need to be tested against experimental data.
Here, we develop and apply these approaches to analyse trap morphogenesis in U. gibba. We show that after forming a near-spherical shape, U. gibba traps undergo defined changes in shape and curvature. By measuring 3D snapshots of traps at various developmental stages and exploring computational growth models, we show that differential rates and orientations of growth are both likely involved in the observed shape transformations. This hypothesis is further tested by marking cells with green fluorescent protein (GFP) and testing the resulting cell and clone shapes against model predictions. To account for oriented growth, the computational model invokes a proximodistal polarity field that is comparable to that proposed to account for morphogenesis during these stages involves transformation of a near-spherical overall shape to an oblate spheroid (a sphere that has been squashed along one axis) with a straightened edge.
At very early phases of development, traps were hidden from view because they were held within a spiral structure, termed the circinate apex [39]. We imaged traps at daily intervals from when they emerged from the spiral until they reached maturity. Trap length was estimated according to the distance from the dorsal lip ( Fig 1C, red square) to the furthest point at the back of the trap (Fig 1A, red line).
Plotting log of trap length against time gave an estimated strain rate (relative growth rate) of 1.8% h −1 ± 0.13 ( Fig 3J, S3 Data, all measured rate estimates are given with ± twice the standard deviation). The growth curve was extrapolated back in time (Fig 3J, dotted line) to define an initiation time (0 days after initiation [DAI]), corresponding to a length of 10 μm (i.e., approximately 1-2 cells). Using this growth curve, a standard time in DAI could be assigned to any trap based on its length (blue region, Fig 3J, S4 Data).
The above framework allowed us to determine strain rates for various trap domains that could later be used to constrain parameters in growth models. We first measured circumferences in the three section planes at different stages of development (Fig 3K, S3 Data). The strain rate was higher along the sagittal circumference (1 .65% h −1 ± 0.06) compared to the other circumferences (1.52% h −1 ± 0.07 and 1.39% h −1 ± 0.06). The sagittal section was further divided into three subdomains based on three landmarks that could be identified throughout development (Fig 1C and Fig 3A-3I, green, magenta, and red squares). These landmarks allowed three domains to be defined: ventral midline, dorsal midline, and stalk diameter ( Fig  1C, magenta, red, and green domains). Strain rates for these regions were then estimated from the staged traps ( Fig 3L, S3 Data). The ventral midline grew faster (2.07% h −1 ± 0.09) than the dorsal midline (1.71% h −1 ± 0.09), and the stalk diameter grew the slowest (0.77% h −1 ± 0.13).

Tissue-level modelling
To explore hypotheses that might underlie the observed morphogenetic changes, we developed a series of models constrained by the experimental growth rate data. Models were kept as simple as possible, with hypothetical factors (for example, Midline factor [MID], Stalk factor [STK], Ventral factor [VEN]) being successively introduced to give a clear indication of what contributes to the resulting shape changes. Although model parameters were constrained by experimental data, changes in tissue curvature generated by the model were not specified but were an emergent property arising through mechanical constraints of tissue connectivity.
We used the growing polarised tissue (GPT) modelling framework, in which tissue is treated as a continuous sheet of material with defined thickness, termed the canvas [40]. It is assumed that for each region of the tissue, there is a specified rate of growth that defines how much that region would grow in mechanical isolation from neighbouring tissue. This rate of specified growth is a tensor quantity representing the possibility that the growth (strain rate) may be by different amounts in different directions. Resultant growth is how each region grows in the context of mechanical constraints arising from connectivity with other regions and includes anisotropies, rotations, and curvature that emerge from such constraints [41]. Specified growth, therefore, refers to the intrinsic or active properties of a region, which may be influenced by local gene expression, while resultant growth also includes the passive changes that arise through mechanical connectivity with other regions. Either type of growth can be isotropic (equal in all directions) or anisotropic (greater in some orientations than others). Regional factors in the canvas can modulate specified growth rates, allowing various patterns of growth to be established.
Computationally, the problem is to calculate the deformation field or resultant growth (i.e., a mapping of each region of the tissue to its new position) that will result from applying the field of specified growth rate for all regions when mechanically connected together over some small time interval. In general, there will be no deformation field in which every region of the tissue achieves its specified growth. The difference between the specified and resultant growth is the residual strain and produces a proportionate residual stress. The actual deformation resulting from the field of specified growth is taken to be whatever shape minimizes the residual strain energy. Residual strain is assumed to dissipate after each growth step, reflecting the irreversible plastic flow involved in plant growth [42].
Changes in curvature can arise if regions of the tissue are specified to grow at different rates and/or directions. For example, if a region with high specified areal growth rate is surrounded by regions with low specified growth rate, the tissue may buckle, with the faster-growing region bulging out. More recently, another type of buckling has been described in which specified areal growth rates can be uniform but in which directions of specified anisotropic growth vary [43]. To distinguish between these two causes of buckling, the terms areal and directional conflict resolution have been proposed [44]. In both cases, changes in curvature arise through differential growth properties between regions, with conflicting growth patterns preventing each region from attaining its specified growth rate. Buckling or curvature helps reduce this conflict (reduces strain energy) by allowing each region to grow more closely to its specified rate. Areal conflicts arise when there are nonuniform specified areal growth rates (which may be isotropic or anisotropic), whereas directional conflict can arise when there are nonuniform orientations of growth (even if areal growth rates are uniform). Directional conflicts, therefore, necessarily require specified anisotropy, whereas areal conflicts do not. We explore the ability of each type of tissue conflict resolution, areal or directional, to account for the observed growth rates and shape changes of the U. gibba trap.
In all models that follow, the initial canvas was a hollow sphere with a uniform wall thickness of 30 μm and a diameter of 100 μm, corresponding to the approximate trap shape at 4 DAI ( Fig 2K-2T and Fig 4A). To assist with visualisation in three dimensions, a grid of latitudes and longitudes was superimposed on the initial canvas. All regional factors, as well as the gridlines, were fixed to the canvas and deformed with it during growth.

Shape transformation through areal conflict resolution
We first explored models in which specified growth rates in the plane of the canvas are equal in all directions (isotropic) but can differ between regions, creating potential areal conflicts. Starting from an initial near-spherical shape, one of the key changes during trap morphogenesis is formation of elliptical circumferences in the transverse and frontal planes (Fig 2H and  2I). To a first approximation, this represents a transformation from a sphere to an oblate spheroid. In an oblate spheroid, the two elliptical circumferences are shorter than the circular circumference. In principle, transformation from sphere to oblate spheroid could arise if specified areal growth rate is faster along one circumference, causing the sphere to flatten. This faster-growing circumference could correspond to the sagittal circumference of the Utricularia trap because this ends up being longer than the other circumferences and grows at a faster rate (Fig 3K, S3 Data).
To explore this idea, we introduced a factor, MID, expressed along the sagittal circumference or midline of the canvas (Fig 4A and 4B, red). The concentration of MID was set to be highest at the midline and gradually declined away from it. We next needed to constrain parameters in the model according to observed circumferential growth rates. The areal strain rate for a region of planar tissue is the sum of the two linear strain rates in orthogonal directions within the plane. We therefore set the basic specified areal strain rate of our canvas to the sum of the transverse and frontal circumferential strain rates (2.9% h −1 ) according to parameter b planar (Fig 4C; see model description in Materials and Methods for full details). This strain rate was promoted by MID to give a value of twice the sagittal circumferential strain rate (3.3% h −1 ) along the midline, according to parameter p mid . Specified growth rate in thickness was assumed to be uniform throughout the canvas and set to an experimentally determined average of 0.5% h −1 , according to parameter b thickness (S2 Fig, S10 Data). Thus, the model had a total of three parameters (b planar , p mid , b thickness ) controlling specified growth rates, each of which was experimentally constrained. These constraints were sufficient for the model because specified growth rate was assumed to be uniformly affected by a given factor.
Running this model led to a transformation of the initial sphere to an oblate spheroid ( Fig  4E, S3 Movie). The shape broadly matched that of the front and top views of the mature trap but lacked the straight ventral edge seen in side view (Fig 4D, compare to sagittal shape outline in red). The resultant areal strain rate in the central MID region was slightly lower than that specified (3.18% h −1 compared to 3.3% h −1 ). This deficit arises because the slower growth of the rest of the canvas constrained growth of the MID region, even with the change in curvature (i.e., the areal conflict was not fully resolved). For this reason, in subsequent models, we adjusted parameter values for specified growth rates by trial and error in such a way that they gave resultant strain rates that matched experimental strain rate measurements when the model was run. Unresolved areal conflict also introduced slight resultant anisotropy in the pattern of growth, indicated by the field of maximal growth orientations ( Fig 4F). These differences between specified and resultant growth highlight emergent features arising through mechanical constraints.
In the above model, the region of the midline that intersected with the stalk was specified to grow at the same rate as the rest of the midline. However, the strain rate for stalk diameter was measured to be lower than other midline regions (0.77% h −1 ; Fig 3L, S3 Data). We therefore introduced an additional factor, STK, at the 'South Pole' of the canvas (Fig 4G, green), which inhibited specified areal strain rate according to parameter h stk (Fig 4H). The result of running this four-parameter model (b planar , p mid , b thickness , h stk ) was an oblate spheroid with a slight inflexion at the STK domain caused by the areal conflict between the slower-growing STK domain and its surroundings (Fig 4I and 4J; compare to mature trap outlines shown in red, orange, and cyan and S4 Movie).
In the output of the above model, the mouth region remained close to the stalk, in contrast to the observed displacement of the mouth at later stages (Figs 2G and 3I). This displacement of the mouth reflects the higher strain rate of the ventral midline of 2% h −1 (Fig 3L, S3 Data.). To account for these observations, we introduced a ventral midline factor, VEN (Fig 4K,  magenta), which promoted specified areal strain rate according to parameter p ven (Fig 4L), giving a resultant strain rate along the ventral midline of approximately 2% h −1 . The result of running this five-parameter model (interactions summarised in Fig 5A) is a trap with a longer ventral midline that bulges out (Fig 4M and 4N [red] and S5 Movie), unlike the straightened ventral midline of the real trap (Figs 4D and 1). Bulging arises through the areal conflict caused by the ventral midline region growing faster in all directions than its surroundings. Thus, an areal conflict model can account for the overall flattening of the sphere to create an oblate spheroid but does not readily account for shape of the ventral midline.

Shape transformation through directional conflict resolution
To determine whether a model based purely on directional conflict resolution might better account for the observed morphogenetic changes, we kept specified areal growth rates uniform (3% h −1 , defined by a slightly modified value of b planar ) and varied specified rates of growth in different orientations (specified anisotropy). Although specified areal growth rates were uniform, resultant areal growth rates need not have been because directional conflicts can lead to some regions growing faster or slower than the rate specified. Specified growth rate in thickness (b thickness ) was the same as for the areal conflict model.
To achieve specified anisotropy, we incorporated a polarity field into the model by introducing a source (plus-organiser [+ORG]) and sink (minus-organiser [−ORG]) of a propagating factor, polariser (POL). Taking the local gradient of POL allowed a polarity field to be specified. In a previous model of planar leaf development of Arabidopsis, the +ORG was positioned at the base of the leaf primordium, leading to a proximodistal gradient in POL [45]. By analogy, we introduced a +ORG in the stalk region of the Utricularia trap by having production of POL promoted by STK (Fig 6A). To provide a distal/marginal anchor for the polarity field, we also introduced a −ORG in the mouth region through enhanced degradation of POL (Fig 6A). Following a period of diffusion, POL concentrations were fixed to the initial canvas. Two local specified growth rates (specified strain rates) could then be defined: specified growth rate parallel to the polarity (K par ) and specified growth rate perpendicular to the polarity (K per ).
To generate the transformation from sphere to oblate spheroid, K par was promoted by MID according to parameter p mid (and K per correspondingly reduced to keep specified areal growth rate at 3% h −1 ). The specified anisotropy, (K par − K per )/(K par + K per ), is shown colour-coded in for directional conflict model, with MID and VEN promoting specified growth rate parallel to the polarity and STK inhibiting specified growth rate parallel to the polarity. Note that to maintain constant specified areal growth rate, an increase in K par has to be compensated for by a corresponding decrease in K per (indicated by mutual inhibition). (C) Integrated model. Regulation of K par and K per is separable. MID promotes K par . STK inhibits both K par and K per , and VEN promotes K par and inhibits K per . There is also a further parameter that influences the width of the VEN domain (t ven ) (not shown). (D) KRN for Arabidopsis leaf model for comparison. The MID factor for the Arabidopsis model is expressed in the midline region and has a higher level of expression in the proximal half of the primordium. The LAM factor is expressed in the presumptive lamina, which occupies most of the primordium except for its most proximal region. PGRAD has a graded distribution that decreases from proximal to distal positions. LATE is expressed uniformly and increases with time. https://doi.org/10.1371/journal.pbio.3000427.g005 Shaping of a three-dimensional carnivorous trap through modulation of a planar growth mechanism Fig 6B. This three-parameter model (b planar , p mid , b thickness ) led to an oblate spheroid shape ( Fig  6C, S6 Movie). The oblate spheroid was narrower than in the areal conflict model as the midline region grew less in width and was a better fit when compared to mature trap shape frontal (orange) and transverse (cyan) outlines. Also, in contrast to the areal conflict model, the maximal rate of resultant growth within the midline regions was oriented parallel to the midline rather than perpendicular to it (compare Figs 6D and 4F).
We next inhibited K par by STK according to parameter h stk (Fig 6E and 6F). This fourparameter model (b planar , p mid , b thickness , h stk ) gave an oblate spheroid with a wide STK domain ( Fig 6G and 6H, S7 Movie). The STK domain grew in width because of the high value of K per needed to keep the total specified areal strain rate constant. To generate an elongated ventral midline, K par was further promoted by VEN according to parameter p ven (Fig 6I and 6J). This five-parameter model (interactions summarised in Fig 5B) gave an extended ventral midline region that was relatively straight compared to the bulged-out shape generated by the areal conflict model (compare Fig 6K and 6L, S8 Movie, with Fig 4M and 4N, S5 Movie). Thus, the directional conflict model accounted for the main shape transformations of the trap more effectively than the areal conflict model. However, the directional conflict model gave a less rounded STK domain than the areal conflict model and thus matched this aspect of development less effectively.

Shape transformation through integrated areal and directional conflict resolution
To determine whether the various features could be captured with a single model, we developed an integrated model incorporating both directional and areal conflicts (seven-parameter model; interactions summarised in Fig 5C). To achieve this, we removed the constraint from the pure directional conflict model that specified areal growth rates were uniform. This allowed the specified areal growth rate in the stalk region to be lower than the rest of the canvas and enhanced in the midline, similar to the areal conflict model (compare Fig 6M with Fig   Fig 6. Tissue-level modelling of trap development through directional conflict resolution. Specified areal growth rate is uniform in A-L (directional conflict model) but not in M-P (integrated model). Initial spherical canvas is shown from side and front in the left two columns. Resultant shapes from side, front, and top are shown in right three columns, with experimentally observed circumferences (red, orange, cyan as shown in Fig 4D) Fig 4M). (L) Midsection through resultant shape (contrast with Fig 4N). Shaping of a three-dimensional carnivorous trap through modulation of a planar growth mechanism 4L). We also broadened the domain of anisotropy in the VEN region (compare Fig 6N with  Fig 6J), allowing for a better match to the cell-shape data described below (Fig 6P). These changes involved introducing two further parameters controlling specified growth rates. Running this model gave a shape showing a good match to that observed, with a small rounded STK domain at the final stage ( Fig 6O and 6P, S9 Movie). A similar shape was generated if variation in rate of growth in trap wall thickness was incorporated, showing this did not have a marked effect (S3 Fig).Thus, the integrated model could account more effectively for observed trap shape transformations than models based purely on directional or areal conflict resolution (compare Fig 6O and 6P, S9 Movie, Fig 6K and 6L, S8 Movie with Fig 4M and 4N, S5 Movie).
Comparing the Growth regulatory network (KRN) of the integrated model ( Fig 5C) with that proposed for Arabidopsis leaf development (Fig 5D, [45]) reveals both similarities and differences. In both cases, factors expressed in the midline region inhibit K per . In Utricularia, the factor is VEN, whereas in Arabidopsis, it is MID, which is expressed most strongly in the proximal midline. In both cases, K per is low in the basal part, leading to a narrow supporting structure (stalk or petiole). In Utricularia, this is implemented by STK inhibiting K per , whereas in Arabidopsis, it is through a lamina factor (LAM) promoting K per . In Utricularia, STK also inhibits K par because STK defines a domain intersecting the stalk and thus affects both its width and thickness.
A notable difference between Utricularia and Arabidopsis is that K par is promoted in the midline regions of Utricularia (by VEN and MID) but not in Arabidopsis. This difference reflects the planar nature of Arabidopsis leaf growth. If the Arabidopsis midline region grew faster in length than the adjacent lamina, the midline would buckle out of the plane. In Utricularia, in which planarity is not required, enhanced growth of the midline regions leads to the oblate spheroid shape and increased length of the ventral midline. A further difference between the species is modulation of growth by a temporally varying factor (LATE) and a graded proximodistal factor (PGRAD) in Arabidopsis. In the absence of live imaging of regional growth in Utricularia, it is not clear whether such factors may also be involved in Utricularia trap development.
A further feature of the U. gibba model is that it allows evolutionary variation in trap shape to be explored. Utricularia traps vary in shape between species from terminal types, which have the mouth distant from the stalk (for example, U. bisquamata, Fig 7A), to basal types that have the mouth positioned near the stalk (for example, U. praelonga, Fig 7B) [46]. U. gibba belongs to a lateral type, intermediate between these extremes ( Fig 7C). To illustrate how the model could be modulated to generate these forms, we varied the parameters by which MID and VEN affect K par . Increasing promotion of K par by VEN and reducing promotion by MID gave a shape resembling the terminal type ( Fig 7D), whereas decreasing promotion of K par by VEN and increasing promotion by MID gave the basal type ( Fig 7E). The resultant shape of the integrated model is closest to the U. gibba intermediate type ( Fig 7F). The extent to which such variations are valid could be tested by analysing the growth of each trap type.

Cellular-level data and modelling
To explore how cell growth and division may be integrated within these models, we tiled the canvas with virtual cells (polygons) with vertices that are displaced as the canvas grows [6]. New walls may then be introduced through cell division as cells reach a threshold size. In this modelling framework, cell divisions do not contribute to specified growth rate of the tissue but respond to the pattern of local growth. This is consistent with the mechanism of plant cell growth, which is driven by turgor pressure stretching the cell walls. Cell divisions thus provide new partitions that maintain the material properties of the growing cell-wall mesh rather than being drivers of growth.
To provide constraints for such cellular-level models, we first estimated cell division parameters by counting cell numbers along different circumferences and regions (Fig 8A-8F). For all regions except the stalk, cell numbers showed an approximately exponential increase until about 6-7 DAI, after which cell numbers levelled off ( Fig 8G-8L, S5 Data). This suggests an early phase when cell divisions occur, followed by a later phase after 6-7 DAI, when cell division slows down or arrests. During the division phase, the rate of increase in cell number was 1%-2% h −1 , comparable to the strain rates ( Fig 3K and 3L, S3 Data), suggesting cell division broadly keeps up with growth. The reduced rate of division after 6-7 DAI did not correlate with a change in growth rate, which reduced later (approximately 11 DAI, Fig 3K and 3L, S3 Data), consistent with division not being the driver of growth. Estimates of cell area during the division phase gave a mean of approximately 50 μm 2 (S4 Fig, S6 Data). The stalk region showed very little change in cell number, indicating that cell division rates were low in this region over the period analysed ( Fig 8L, S5 Data). This lack of division correlated with a slow growth rate.
Cellular-level models were developed according the principles established in the analysis of Arabidopsis leaf development [6] but extended from a flat sheet to a curved sheet embedded in three dimensions. The surface of the initial canvas was tiled with an array of cells with mean cell area of 50 μm 2 (approximately 350 cells). As in Arabidopsis, we invoked a dual-control model in which there is spatiotemporal regulation of both growth (K par and K per ) and cell division. Cell division required expression of a factor conferring division competence (CDIV) that was expressed throughout the canvas (except in the STK region) until 6.5 DAI, after which it was switched off. Execution of division occurred when cells reached a threshold area, T A , set to  Assuming symmetric division, this would give daughters of 35 μm 2 and a mean cell area of approximately 50 μm 2 . Thus, incorporating cells involved two additional parameters, both of which were experimentally constrained: mean division threshold T A = 70 μm 2 and inactivation of CDIV in all cells at 6.5 DAI. For simplicity, the mouth region was modelled in a similar way to the rest of the tissue and is shown in grey in output images. Cell sizes and cell shapes were an emergent property that resulted from running the model.
Running the areal, directional and integrated models with these assumptions gave cell areas in the main trap of less than 70 μm 2 until 6.5 DAI, after which cell size increased following division arrest (integrated model shown Fig 9A, S10 Movie, S11 Movie [other models in S5 Fig, S12 Movie and S13 Movie]). Cells in the STK region enlarged from an earlier stage because they did not divide because CDIV was absent.
To evaluate these models against experimental data, we segmented 3D confocal images of traps at corresponding stages to those shown for the model (Fig 9B). The observed cell sizes showed more variation in spatial pattern than generated by our simplified models, indicating that more elaborate mechanisms operate for spatiotemporal control of division and/or growth than those implemented. Nevertheless, broad trends could be compared. As with the model outputs, cell area increased after 6-7 DAI (S4 Fig, S6 Data), except for a subgroup of cells (hemispherical gland cells) that remained small (S6 Fig). Final cell areas in the experimental data showed no major enhancement or reduction in cell areas in the ventral midline (arrowed in the two examples of Fig 9C and 9D). By contrast, in the areal conflict model, cells were larger in the midline regions because these were the regions of higher growth rate (arrowed in Fig 9E). In the directional conflict model, cells along the ventral midline were smaller (arrowed in Fig 9F). This reduction in cell size arises through the directional conflict, which leads to reduced resultant areal growth rate in this region throughout the simulation. This effect was absent in the integrated model because K per was reduced by less and over a broader domain (arrowed in Fig 9G). Thus, the integrated model gave the best overall match to the pattern of cell sizes in this region (compare Fig 9G to Fig 9C and Fig 9D).
Another possible test of the models might be the pattern of cell-shape anisotropy because the models make very different predictions. The areal conflict model predicts resultant growth anisotropy perpendicular to the midline (Fig 4F), whereas the directional conflict model predicts resultant anisotropy parallel to the midline, particularly in the ventral midline ( Fig 6D). These anisotropies in growth will affect cell shape after cell divisions cease. To evaluate these effects, we colour-coded the cells generated by models according to their cell-shape anisotropy (defined as R − 1/R + 1, where R is the ratio of the long/short axis of an ellipsoid fitted to the cell) and showed the orientation of the cell long axis with a black line for cells with strong anisotropy (results for integrated model shown in Fig 10A, other models in S7 Fig). To allow direct comparison with the experimental data, we used the same colour-coding system for the segmented cells of traps ( Fig 10B).
The most striking region of high anisotropy in the experimental data was in the ventral region, where the long axis of the cells was oriented mainly parallel to the ventral midline (arrowed in the two examples showing in Fig 10C and 10D). This finding was inconsistent with output of the areal conflict model (Fig 10E, S14 Movie) but was predicted by the directional conflict model (Fig 10F, S15 Movie). An even better match with the experimental data Shaping of a three-dimensional carnivorous trap through modulation of a planar growth mechanism was obtained with the integrated model ( Fig 10G, S16 Movie, S17 Movie), which invoked a broader domain of specified anisotropy in the ventral midline.
In principle, further modifications and parameters could be incorporated to the integrated model to give a better match to the distribution of cell patterns or shape of the trap. Additional features, such as introducing a discontinuity at the mouth and simulating its opening at later stages, could also be introduced. However, without further experimental data to test and constrain the modelling, such an exercise may not be mechanistically informative.

Clonal analysis
Cell-shape analysis provides evidence for growth anisotropy after cell divisions have ceased (after 6-7 DAI). To determine whether growth anisotropy was also present during earlier Shaping of a three-dimensional carnivorous trap through modulation of a planar growth mechanism phases, when divisions were occurring, we used clonal analysis to examine the size and shape of patches of cells generated by single progenitor cells. We achieved this by developing a transformation protocol for U. gibba (detailed in Materials and Methods) and introducing a heatshock (HS)-inducible Cre-lox system generating clones expressing GFP on an mCherry background.
We first analysed the pattern of the shapes of clonal regions generated by each model. Virtual clones were generated by colour-coding approximately 40% of cells at approximately 4 DAI and following their descendants (Fig 11). Virtual clones could have two components contributing to their overall anisotropy in shape. First, there could be more cells along the major axis of the clone (cell-number anisotropy), arising because of anisotropic growth during the Shaping of a three-dimensional carnivorous trap through modulation of a planar growth mechanism period when cell divisions occur. Second, cells could be elongated along the major axis of growth (cell-shape anisotropy), arising because of anisotropic growth after division arrest. The anisotropy of the clone as whole (clone anisotropy) is the product of cell-number and cellshape anisotropy. In the directional and integrated models, both types of anisotropy contribute to the elongated clones of the ventral midline region (Fig 11F and 11J arrowed, S18 Movie, S19 Movie [as compared to the areal model, S20 Movie]).
To compare these predicted patterns to clones generated experimentally, we induced clones using the Cre-lox system [47]. We introduced a construct with the cauliflower mosaic virus (CaMV) 35S promoter driving GFP interrupted by an mCherry coding sequence with a terminator, flanked by lox recombination sites (see Materials and Methods). The construct also carried Cre recombinase under the control of an HS promoter. Following HS, GFP sectors were visualised 4 days later, when the trap size was at that expected for approximately 10 DAI. Clones between 3-30 cells were selected for measurement. Fig 11M-11P (and S21 Movie, S22 Movie) illustrate some of the sectors obtained. A summary of results for clones in the different regions of the trap are shown in Fig 11Q-11T, (S17 Data). In the ventral midline region, clones were preferentially oriented parallel rather than perpendicular to the midline, consistent with the directional conflict and integrated models.
To quantify these components of anisotropy in the experimental data, we first subdivided the trap into ventral midline, dorsal midline, and lamina domains. A clone was considered to be within the midline domain (dorsal or ventral) if its centre was a distance of five cells or less from the midline. For each domain, we measured the ratio of the long/short axis of each clone (clone anisotropy), the ratio of cell number along the long axis and short axis of each clone (cell-number anisotropy), and clone anisotropy divided by cell-number anisotropy (cell-shape anisotropy). Clone anisotropy was significantly higher in the ventral midline region compared to the other regions ( Fig 11U, S7 Data). Some of this difference came from cell-number anisotropy, which was significantly higher in the ventral compared to dorsal midline regions ( Fig  11V, S7 Data). Cell-shape anisotropy was also significantly higher in the ventral midline region (Fig 11W, S7 Data). These results are more consistent with the directional conflict and integrated models than the areal conflict model (compare Fig 11E-11H with Fig 11I-11L and Fig  11A-11D and S18 Movie, S19 Movie, and S20 Movie) and indicate that anisotropic growth occurred in the ventral midline region both during the phase of cell division and after cell division arrest.

Evidence for a polarity field
The above models involving directional conflicts hypothesised a polarity field running from stalk to mouth. To test whether such a polarity field exists, we analysed the pattern of glands on the inside of the trap because hair morphology has been used to infer cell polarity fields in several cases. As previously noted [13], quadrifid glands, which decorate the inside of the mature trap, often have arms more splayed out at one end than the other (Fig 12, S23 Movie).
To determine whether this polarity of the quadrifids is organised as a field, the glands were imaged in three dimensions using OPT and confocal microscopy.  Shaping of a three-dimensional carnivorous trap through modulation of a planar growth mechanism were identified on each quadrifid, one at the centre and one at each arm tip (Fig 12A-12C). The landmarks were placed by examining each quadrifid in isolation. The two pairs of nearest tip landmarks were identified and the distances between them calculated and subtracted from one another. The quadrifid was assigned a polarity, pointing towards the more widely splayed arms and passing through the middle of the quadrifid (Fig 12D, arrow). A threshold value was applied to identify glands in which polarity could not be assigned with confidence. In these cases, quadrifid was assigned an axiality (line passing through the middle of the quadrifid) but no polarity. Quadrifid polarities and axialities were viewed in three dimensions (S24 Movie, S25 Movie) or projected in 2D (Fig 12). An additional measure was applied to confirm that quadrifid polarities based on distance between arms was consistent with relative arm lengths (S8 Fig).
In side views, quadrifid polarities mainly pointed away from the stalk region and towards the mouth (for example, 35/37 in the example shown in Fig 12E and 12F, S24 Movie, S8 Data, another example shown in S8D Fig), consistent with the hypothesised polarity field (Fig 12G). Ventral (Fig 12H and 12I, S25 Movie, S8G Fig, S8J Fig) and dorsal views (Fig 12K and 12L,  S8M Fig, S8P Fig) also indicated polarity running from stalk to mouth, consistent with the proposed polarity field (Fig 12J and 12M). The sense in which the polarity arrow points is arbitrary and could equally well be depicted as going from mouth to stalk. The key observation is that the polarity is coordinated to point preferentially in one direction over the other. Thus, the quadrifid analysis supported the polarity field hypothesised in the model.

Discussion
We show that following formation of a near-spherical shape at early stages of development, the Utricularia trap flattens to form an oblate spheroid and develops an extended ventral midline. We further show that this shape change, as well as the broad pattern of cell shapes and sizes, can be accounted for by a model that invokes many of the same principles as those used to account for Arabidopsis leaf development but in the context of a highly curved tissue sheet. Our results thus indicate that simple modulations of core developmental processes can account for shape changes in a diverse range of contexts.
The model we propose has experimentally constrained parameters influencing specified growth rates (parallel or perpendicular to a polarity field) and affecting cell division (threshold area for division execution and duration of division competence). Changes in organ curvature and cell sizes and shapes are not directly specified but emerge from running the model because of mechanical connectivity between regions. By introducing different elements into the model progressively and comparing outputs to experimental observations, the contribution of particular components could be evaluated.

Organisation of polarity field
To account for Utricularia trap morphogenesis, we invoked a proximodistal polarity field to provide orientational information for specified anisotropic growth (Fig 6). The polarity field radiates out from a +ORG at the base (stalk) and converges distally towards a −ORG near the mouth region. The field is comparable to that proposed for Arabidopsis leaf development, Shaping of a three-dimensional carnivorous trap through modulation of a planar growth mechanism except that instead of propagating on a flat sheet, it propagates through a curved sheet. The polarity field is parallel to the midline in both cases. However, in contrast to the Arabidopsis leaf, in which the midline is a linear continuation of the petiole, the midline of the Utricularia trap forms a curve that traverses the stalk, reflecting its peltate organisation [16].
Evidence for the proposed polarity field in Utricularia comes from the analysis of quadrifid glands. These glands can be used as indicators of polarity, similar to hairs in Drosophila wing [48] or trichomes on leaves of Arabidopsis or barley [49][50][51]. The coordinated polarity exhibited by these hairs, glands, or trichomes likely reflects polarity fields established much earlier in development, for example, by planar or tissue cell polarity pathways [43,[52][53][54][55][56][57].
In our model, we propose that a key role of the polarity field in Utricularia is to coordinate the orientation of specified anisotropic growth. Alternatively, specified anisotropic growth could be coordinated through the orientation of stresses [58]. Stresses can originate from differential growth of a tissue sheet (areal conflict), and it has been proposed that such residual stresses might change local cell-wall properties to orient specified anisotropic growth without the need for a polarity field [59]. However, computer simulations show that such a model does not allow local orientations to be specified in a stable manner [41,59]. Another possibility is that differential pressure within the trap generates circumferential stresses that orient growth. In normal conditions, the interior of the mature trap is under low pressure caused by water being pumped out of the trap chamber [36,60,61]. The differential pressure could create circumferential stresses that orient growth. However, under the tissue culture conditions employed for the growth of Utricularia here, the traps develop with little pressure differential (their shape is similar to the relaxed state), making this mechanism unlikely.

Contributions of directional and areal conflict resolution
We show that both areal and directional conflict resolution play a role in shaping the Utricularia trap. A role for areal conflict resolution is suggested by the slow growth rate in diameter of the stalk. A role for directional conflict is suggested by tissue-level modelling: generation of an extended ventral midline is most readily achieved through greater specified growth parallel compared to perpendicular to the midline (i.e., anisotropic-specified growth). We found that with isotropic-specified growth alone, the ventral midline bulges out. Further support for specified anisotropy derives from clonal analysis, which indicates greater cell proliferation parallel compared to perpendicular to the ventral midline. Cell-shape analysis also shows that cells are more elongated parallel to the ventral midline, consistent with anisotropic growth after the cessation of division. Enhanced growth parallel to polarity along the midvalve region of Capsella fruits (which corresponds to the midline of the carpel primordia) has also been proposed to account for shape change [62].
Comparisons between the growth model of the Utricularia trap and that proposed for Arabidopsis leaf development reveal both similarities and differences. In both cases, specified growth perpendicular to the polarity is inhibited in the midline and proximal domain, leading to a narrow midline and supporting structure (stalk or petiole). However, growth parallel to the polarity is promoted in the midline regions of Utricularia but not in Arabidopsis. This difference reflects the planar nature of Arabidopsis leaf growth. If the Arabidopsis midline region grew faster in length than the adjacent lamina, the midline would buckle out of the plane. In Utricularia, where planarity is not required, enhanced growth of the midline regions leads to the oblate spheroid shape and increased length of the ventral midline.

Cell division and growth
During the early phases of U. gibba trap morphogenesis, cell division occurs concurrently with growth, whereas at later stages, growth occurs in the absence of division, leading to cell expansion. This is comparable to the situation in planar leaf development [63]. We model cell division in the U. gibba trap through a dual-control mechanism in which both cell division and growth parameters are under spatiotemporal regulation. We assume cells are competent to divide at early stages (before 6.5 DAI), after which divisions are arrested. Competent cells execute division when they reach a threshold size. Such a model broadly accounts for observed clonal anisotropy in both cell number and cell shape in clones. Anisotropy in cell number arises through anisotropic growth during the cell division phase, whereas anisotropy in cell shape arises through anisotropic growth after divisions have ceased. This type of model contrasts with that proposed for Sarracenia pitcher leaf development, in which changes in the planes of cell division are proposed to drive the development of leaf shape [64]. In our model, the planes of cell division follow as a result of differentially oriented growth rather than being the primary cause of morphogenesis.

Limitations of the model
The model we present aims to capture overall external shape change of the trap but does not account for the more subtle patterns of spatiotemporal variation in cell sizes and shapes of the trap. It also does not account for formation of the mouth opening or the internal foldings and thickenings that occur within the trap, generating the threshold and trap door (Fig 1B and  1C). The model considers only the later stages of morphogenesis after the trap is a near-spherical shape. Curvature of the trap primordium is present before this stage, appearing as soon as the sheet-like nature of the primordium becomes evident. This suggests that sheet formation, which depends on control ab/adaxial patterning in planar leaves [4], is intimately linked with formation of the spherical shape. Further hypothesis development, including the modelling of sheet formation [65] and experimental testing, may be feasible in U. gibba, given its small genome and potential as a model system for carnivorous plants [66].

Evolution of leaf shape
The shape of Utricularia traps is highly constrained by the need to have a sealed trap door that allows low pressure to be established within the trap lumen and released upon triggering [27,36,37,67]. Despite this constraint, Utricularia traps vary in shape between species from terminal types, which have the mouth distant from the stalk, to basal types that have the mouth positioned near the stalk [46]. U. gibba belongs to a lateral type, intermediate between these extremes. We show that a simple mechanism for generating these trap morphologies is to vary the relative rates of growth along the ventral and dorsal midlines, with basal types having the lowest and terminal types the highest ratios of ventral to dorsal growth. If specified growth was purely isotropic, changing these ratios may have consequences on trap function by modifying trap curvature around the midlines. Our analysis shows that orienting anisotropic growth through a polarity field may allow such changes in curvature to be controlled by varying the extent of growth parallel and perpendicular to the polarity. Thus, the morphogenetic system described here may provide great flexibility in allowing shape to evolve even when under strong functional constraints.

Plant material and specimen preparation
U. gibba (The Fly Trap Plants, Bergh Apton, UK) was grown in a heated glasshouse set to 22˚C in plastic trays with 3 cm 1/1 peat/silver-sand mix, topped up with 500 ml deionised water. After flowering, collected seeds were surface sterilised for in vitro culture in 70% ethanol, 0.1% SDS solution for 5 minutes; washed in water, 1 × 4% parazone bleach, 0.2% triton X-100 (11332481001, Merck, Darmstadt, Germany) treatment for 10 minutes; and washed three times in sterile water. Seeds were germinated in sterile pots (100 ml Sterilin jar, 185AM; Slaughter, Basildon, UK) with 25 ml of solid 1B culture medium (

Growth tracking
A cut piece of U. gibba stolon approximately 3 cm long from in vitro culture was placed in liquid 1B medium in a small Petri dish (Sterlin, 50-mm diameter, 124; Slaughter). The youngest trap after emergence from the circinate apex (approximately 150 μm) was imaged every 24 hours until maturity under bright-field light on a Leica M205C stereomicroscope with a Leica DFC495 camera (Leica, Milton Keynes, UK). Trap length (Fig 1A) was measured (Leica LAS version, 4.2 software) and natural log of length plotted against time (Microsoft Excel with LINEST function, Fig 3J, S3 Data). The growth curve was extrapolated back to when the trap was 10 μm long, corresponding to 1-2 cells, which we took to be the initiation stage.

Imaging trap morphology and quadrifid glands
Propidium iodide (PI) staining was applied to fixed U. gibba traps and circinate apexes to achieve maximum depth visualisation [68]. For confocal microscopy, tissue was mounted on cavity slides with an additional gasket for added depth (Frame-Seal Incubation Chambers, SLF0601; Bio-Rad, Hercules, CA, USA). Whole traps were imaged in sagittal orientation on either a Leica SP5-II or Zeiss LSM780 confocal microscope (excitation 514-nm laser line, detection at 580 to 660 nm; Leica). PI staining was most effective in traps <400 μm in length, above which tissue damage occurred. It was possible to get full 3D scans of PI-stained traps up to 115 μm in length. Traps larger than this were imaged to half-trap depth, ending at the dorsal midline vein.
OPT was used to visualise the full 3D shape of traps above 200 μm in length. After PI staining, traps were washed in water, embedded in 1% low melting point agarose (UltraPure LMP agarose, 16520; Invitrogen, Carlsbad, CA, USA), dehydrated overnight in methanol, and cleared in 1 part benzyl alcohol (402834; Merck):2 parts benzyl benzoate (B6630; Merck) (BABB) and prepared for OPT with a prototype scanner, as previously described [38]. UV light was used to view PI fluorescence through the Texas red (TXR) exciter filter 560/40 nm, barrier filter 610 LP. White light was used for a transmission OPT channel.
Live traps from transgenic plants containing fluorescence markers (see constructs section) were mounted in water in 1.2-mm cavity slides (BR475505; Merck) and imaged in sagittal orientation with a Leica SP5-II confocal microscope (10× or 20× lenses). GFP was excited at 488 nm and detected between 500-530 nm. mCherry was excited at 561 nm and detected between 575 and 630 nm. To visualise quadrifid glands, live traps containing fluorescence markers were cut in transverse plain with a razor, mounted in water in cavity slides, and imaged with a Leica SP5-II confocal or Zeiss LSM 5 Exciter (Zeiss, Cambourne, UK) microscope.
To view quadrifid gland pattern in three dimensions with OPT, traps were stained overnight in 10 ml water with 20 μl of 2.5% w/v Toluidine Blue (198161; Merck). Stained traps were washed in water, embedded in 1% low melting point agarose (UltraPure LMP agarose, 16520; Invitrogen), and OPT-scanned in water, with no fixation or clearing, using white light transmission OPT on a prototype OPT scanner (Lee and colleagues [38]).

U. gibba development
To allocate traps to developmental stages according to time (DAI or Hours After Initiation [HAI]), confocal and OPT data sets were clipped to the centre in the sagittal plane in Vol-Viewer, and length was measured from the dorsal lip landmark (Fig 1C) to the furthest point at the rear of the trap (Fig 1A). This length was used to place traps on the mean growth rate trend line (Fig 3J) with a Microsoft Excel histogram macro (S4 Data). To account for shrinkage on dehydration and clearing (5.78% ± 0.45 [n = 6], S9 Data), trap length of fixed traps was increased by 5.78% and DAI or HAI calculated as above.
Trap circumferences and regions (Fig 3, S3 Data) and trap thickness (S2 Fig, S10 Data) were manually measured with VolViewer, and cells were counted (Fig 8A-8F, S5 Data) on clipped planes: sagittal plane was clipped to dorsal midline vein, frontal plane was clipped to front of stalk in the tallest central trap region, and transverse plane was clipped to the trap centre between mouth threshold and trap door. Where half-traps were imaged, circumference measurements and cell counts for transverse and frontal planes were doubled (Fig 3K-3L and Fig 8G-8L).
Data were plotted as histograms in Microsoft Excel with LINEST function for trend line growth rate and standard deviation calculations (S3 Data, S5 Data).
To view trap shape circumferences in each plane, images from VolViewer were cropped in Adobe Photoshop and ellipses fitted and combined with Adobe Illustrator. For mature traps, OPT images had ellipses fitted in transverse and frontal views (Fig 2E and 2H, Fig 2F and 2I). Shape was traced for the sagittal circumference (Fig 2G and 2J) (mean 899 μm long, 10.4 DAI, n = 6). For young traps, ellipses were applied to confocal trap images (mean 55 μm long, 4 DAI, n = 7) (Fig 2O-2T, S1 Data, and S2 Data).

Impact of triggering on mature trap shape
Mature culture-grown traps imaged using OPT did not exhibit the concave wall shape seen in primed, glasshouse-grown traps in water (Fig 1B, Fig 2E and 2F, S1 Fig). To determine how growth of traps cultured in vitro (in liquid 1B culture medium grown in a controlled environment room at 25˚C, with 16-hour light and 8-hour dark cycles) and staining/clearing for OPT influenced shape, traps were mechanically triggered and imaged. U. gibba traps grown in a glasshouse (22˚C) in water were imaged from above under bright-field light (Leica M205C stereomicroscope with Leica DFC495 camera) before and after mechanical triggering with fine forceps (S1A Fig, S1B Fig). Different live-water-grown traps, with and without triggering, were embedded in 1% low melting point agarose (UltraPure LMP agarose, 16520; Invitrogen), and OPT-imaged in water on a prototype OPT scanner [38] (S1E Fig, S1J Fig). To explore whether dehydration and clearing for OPT caused triggering, traps grown in water and in vitro in liquid 1B medium were first OPT-imaged in water, then dehydrated overnight in methanol and cleared in BABB (Merck) and OPT-imaged again (S1K-S1V Fig). Dehydration and clearing acted to trigger the trap, resulting in rounder shapes (S1O-S1P Fig, S1U-S1V  Fig). Traps grown in water were larger than those grown in liquid 1B medium and showed greater shape change on triggering. To determine whether culture-grown traps could regain the primed shape, traps were imaged under bright-field light on a Leica M205C stereomicroscope with Leica DFC495 camera from above before and 24 h after transfer from liquid 1B culture medium to water, demonstrating recovery of the primed shape of water-grown traps (S1C-S1D Fig, S9 Data).

U. gibba transformation
35S::loxRFPloxGFP-HSP18::CRE-35S::Kan (EC71194) HS-inducible plants were screened for GFP and mCherry fluorescence on a Leica DM6000 fluorescence microscope, and any lines with GFP fluorescence before heat shocking were discarded. Transformed plants were confirmed to be single-copy by iDnaGENETICS, Norwich, UK. Samples were analysed by qPCR using a multiplexed taqMan reaction assaying for NPT2 and the 35S promoter. Cofactor of nitrate reductase and xanthine dehydrogenase (CNX3, Scf00029.g3638.t1) was used as a single-copy control. CNX3 was reported to be single-copy in U. gibba (Ibarra-Laclette and colleagues [22]), and this was confirmed by BLAST analysis within genome assemblies of the Bergh Apton accession used in this study [70].

Clonal analysis
Growing stolon tips of plants, 2-3 cm in length, were collected for HS treatment and placed in six-well plates (657160; Greiner Bio-One LTD, Stonehouse, UK), each well containing 5 ml 1B media (see Plant material and specimen preparation section above) and 4-6 growing tips. Plates were sealed (Micropore tape) and floated in a 45˚C water bath for 6-8 minutes. HS tissue was left to grow under standard in vitro conditions for 4 days. Traps 772-1090 μm in length (10 to 11 DAI) showing GFP clones were selected by visualising them under a fluorescence microscope (Leica Fluo III stereo or Leica DM6000). These traps were imaged in sagittal view and rotated under the coverslip or cut with a razor to allow imaging of GFP clones at multiple angles with a Leica SP5 II confocal microscope (Fig 11M-11P).
Clone anisotropy was calculated by dividing clone length along its longest (major) axis by width along the perpendicular (minor) axis (Fig 11U, S7 Data). Cell-number anisotropy was calculated by manually counting cells along the major and minor clone axes and taking the ratio (VolViewer) (Fig 11V, S7 Data). Cell-shape anisotropy is clone anisotropy/cell-number anisotropy (Fig 11W, S7 Data). To view clone shape and orientation from multiple traps together, clone images from VolViewer were cropped in Adobe Photoshop and ellipses and major axes fitted and placed in their approximate location on a cartoon trap outline with Adobe Illustrator (Fig 11Q and 11T, S17 Data).

Segmentation
Owing to the distinct qualities of images of PI-stained fixed traps and GFP-expressing live traps, two different processing pipelines were applied to extract cellular information from confocal image stacks (Figs 9 and 10).
For GFP-expressing live traps, a Hessian-based membrane enhancing filter [71] was first used to enhance the definition of the outer surface of the trap. Thresholding and morphological operations, followed by a level-set segmentation [72], were used to locate the surfaces of the outermost layer of cells. Triangulated bladder surfaces were extracted from these binary masks using a basic 3D surface net [73], positioning vertices at the voxel centroids. Bilaplacian smoothing was used to give a smoother surface. Considering the signal intensity along a line segment normal to each vertex, vertices were moved to the point with maximum signal intensity in order to better capture the outermost surface of the outer layer of cells. Repeated rounds of bilaplacian smoothing and surface subdivision (splitting each triangular face into four) were used to generate the final refined surface. Stack signal fluorescence was projected onto each vertex of the triangulated surface [74]. Regions occupied by each cell were identified using the Surface Segmentation Potts Model (SurfaceSPM), which extends the method (Segmentation Potts Model [SPM], details to be published elsewhere) to surface image data on triangulated surfaces. The SurfaceSPM is a stochastic procedure, and combining five segmentation runs for each trap with differently-seeded random number generators yielded more accurate segmentations. The SurfaceSPM procedure sometimes generates disconnected cell labels, so labelled regions were divided into connected components, and very small (<0.05 mean label area) components were removed.
Surfaces were clipped using a manually specified polygonal region. Cells touching the edge of the clipped surface or above some size threshold (5 times the mean cell area) removed. For the purposes of cell-number quantification, gland cells were identified as cells with areas smaller than 0.25 times the mean cell area.
Early traps, imaged using PI staining, were hidden within a tight spiral structure. As discussed before, VolViewer was used to identify the region occupied by the trap. Segmentation methods (based on SPM) were used to label the regions occupied by each cell in three dimensions. Labelled regions in the exterior to the trap and within stolons were manually identified and removed.
Through binary erosion, cells protruding from the trap surface were eliminated. Morphological operations and surface nets were used to extract a triangulated surface approximating the outer surface of the trap. This surface was smoothed and translated a small distance inwards along the surface normal and underwent further refinement and smoothing.
Optimized segmentation results, either from the segmentation step of the MARS pipeline [71] (reimplemented by timagetk http://gitlab.inria.fr/mosaic/timagetk) or from the 3D SPM, generated a labelled 3D stack. Triangles of the extracted surface were assigned labels according to the label of the voxel containing their centroid. Following this step, surface label data was processed in the same manner as for the GFP stacks.
Cell areas were calculated as the sum of the areas of the triangles occupied by each cell. Cell anisotropies were calculated using the second moment of area, M, which is a matrix with entries where the integral is over all triangles with the label of the cell and x c is the centroid of the cell. The eigenvalues of this symmetric matrix were calculated, and the anisotropy measure, a, is given by Segmentation software may be found here: https://github.com/jfozard/gibba_analysis.

Modelling framework
All models were produced using the GPT framework [40] with GFtbox software, a MATLAB application from http://cmpdartsvr3.cmp.uea.ac.uk/wiki/BanghamLab/index.php/Software. Models used to generate each figure can be downloaded from http://cmpdartsvr3.cmp.uea. ac.uk/wiki/BanghamLab/index.php/Software or https://doi.org/10.6084/m9.figshare. 8966153.v1, Models.7z archive. Fig 4A-4C, Fig 4E and 4F: Isotropic growth promoted by MID at midline. Prior to growth, factor MID was generated along the midline of the initial canvas and allowed to diffuse from this source with a fixed decay rate. Diffusion was inactivated prior to growth, after which MID concentrations were fixed to the canvas and deformed with it. There are three parameters in the model (b planar , p mid , b thickness ) constrained by linear growth rate measurements. The KRN rate equations are as follows: K par = b planar � pro (p mid , i mid ); K per = K par ; and K nor = b thickness , where b planar = 0.0145 is the basic specified growth rate, constrained in accordance with experimental data to give a specified areal planar growth rate of 0.029 h −1 (2b planar ). p mid = 0.165 is the promotion coefficient of MID on growth such that K per + K par = 0.033 h −1 . b thickness = 0.005 h −1 is the specified growth rate in thickness of the spherical sheet and is set to an experimentally observed average. i mid is level of MID factor at each location in the canvas (established during the initial set up). pro(z, i y ) denotes multiply by (1 + zi y ).

Model descriptions
Because K par = K per , specified growth is isotropic, and there is only areal conflict.   Fig 4G-4J. Fig 6A-6D: Anisotropic growth promoted by MID. As with isotropic model (Fig 4A-4C) there are three parameters (b planar , p mid , b thickness ) constrained by linear growth rate measurements. The KRN rate equations are: K par = b planar � pro(p mid , i mid ); K par = min(2b planar , K par ); K per = 2b planar − K par ; and K nor = b thickness , where b planar = 0.015 is the basic specified growth rate, constrained in accordance with experimental data to give a specified areal planar growth rate of 0.03 h −1 (2b planar ). p mid = 0.35 is the promotion coefficient of MID on growth such that growth rate along the sagittal circumference approximately 0.0165 h −1 in accordance with experimental observations. b thickness = 0.005 as with isotropic models. Note that K par + K per = 2b planar = 0.03 everywhere, so there is no specified areal conflict. Fig 6E-6H: Anisotropic growth as in Fig 6A-6D with growth inhibition by STK. As with isotropic model (Fig 4E-4H), there are four parameters (b planar , p mid , h stk , b thickness ) constrained by linear growth rate measurements. The KRN rate equations are as follows: K par = b planar � pro (p mid , i mid ) � inh(h stk , i stk ); K par = min(2b planar , K par ); K per = 2b planar − K par ; and K nor = b thickness , where h stk = 1.5 is the inhibition coefficient of STK on K par such that the resultant strain rate of the stalk parallel to the midline is approximately 0.0075 h −1 . i stk is level of STK factor. Values of b planar , p mid , b thickness are as in the model for Fig 6A-6D. Fig 6I-6L: Anisotropic growth as in Fig 6E-6H with growth promoted by VEN. As with isotropic model (Fig 4K-4N) there are five parameters (b planar , p mid , h stk , p ven , b thickness ) constrained by linear growth rate measurements. The KRN rate equations are as follows: There are seven parameters in the KRN (b planar , p mid , h stk , p ven , b thickness , h ven , t ven ), the first five of which are constrained by linear growth rate measurements (as in Fig 4K-4N and Fig 5I-5L). The parameters h ven and t ven were adjusted to give a pattern of cell areas and anisotropies similar to those seen in Fig 8C and 8D and Fig 9C and 9D. The KRN rate equations are as follows:  Fig 12E). (A) Clipped OPT sagittal view looking into trap at quadrifid glands on the left-hand wall. Arrows (magenta) orient toward greatest distance between quadrifid arms (DistArms output). Arrowheads were unassigned if distance subtraction value between arm sets was <2 μm (shown as lines). (B) DistArms above threshold (all arrows, green or black); DistArms above threshold and polarity assignment further supported by DistArmsSumArms (black arrows), DistArms below threshold (green lines). 35 Tissue-level modelling of trap development through areal conflict resolution with growth promoted by MID (red) and VEN (magenta) and inhibited by STK (green) (Fig 4K-4N). Scale bar is held constant to show increase in size from initial to resultant shape. Final shape is rotated at the end of the movie to show MID and VEN domains in front view. Note the bulge of the ventral midline region (magenta). MID, Midline factor; STK, Stalk factor; VEN, Ventral factor. (AVI) S6 Movie. Directional conflict resolution model with specified anisotropy promoted by MID. Tissue-level modelling of trap development through directional conflict resolution. Specified anisotropy, defined as (K par − K per )/(K par + K per ) promoted by MID (red) (Fig 6A-6D). Scale bar is held constant to show increase in size from initial to resultant shape. Final shape (oblate spheroid) is rotated at the end of the movie to show MID domain and −ORG (cyan) in front view. The pointed shape at the base of the trap is where +ORG is located. Polarity runs from +ORG to −ORG. MID, Midline factor; −ORG, minus-organiser; +ORG, plusorganiser.
(AVI) S7 Movie. Directional conflict resolution with specified anisotropy promoted by MID and inhibited by STK. Tissue-level modelling of trap development through directional conflict resolution. Specified anisotropy, defined as (K par − K per )/(K par + K per ) is positive in MID domain (red) and negative in STK domain (green) (Fig 6E-6H). Scale bar is held constant to show increase in size from initial to resultant shape. Final shape is rotated at the end of the movie to show MID domain and −ORG (cyan) in front view. Note indented shape at the base of the trap where K par is inhibited by STK. MID, Midline factor; STK, Stalk factor; −ORG, minus-organiser. (AVI) S8 Movie. Directional conflict resolution with specified anisotropy promoted by VEN and MID and inhibited by STK. Tissue-level modelling of trap development through directional conflict resolution. Specified anisotropy, defined as (K par − K per )/(K par + K per ) is positive in MID domain (red), further enhanced by VEN (magenta), and negative in STK domain (green) (Fig 6I-6L). Scale bar is held constant to show increase in size from initial to resultant shape. Final shape is rotated at the end of the movie to show MID domain and −ORG (cyan) in front view. Note the extended ventral midline region (magenta) where specified anisotropy is promoted by VEN. MID, Midline factor; STK, Stalk factor; VEN, Ventral factor; −ORG, minusorganiser. (AVI) S9 Movie. Integrated model movie with specified anisotropy promoted by VEN and MID and inhibited by STK. Tissue-level modelling of trap development through integrated areal and directional conflict resolution (Fig 6M-6P). Scale bar is held constant to show increase in size from initial to resultant shape. Final shape is rotated at the end of the movie to show MID domain and −ORG (cyan) in front view. Note the extended ventral midline region (magenta) where specified anisotropy and growth is promoted by VEN. MID, Midline factor; STK, Stalk factor; VEN, Ventral factor; −ORG, minus-organiser. (AVI) S10 Movie. Integrated model colour-coded for cell area. Cellular-level integrated directional and areal conflict model coloured for cell area. Grey region shows approximate location of mouth ( Fig 9G, S5E Fig, S5F Fig). Scale bar is held constant to show increase in size from initial to resultant shape. Final shape is rotated at the end of the movie. Note similar size of cells in the ventral midline region compared to the rest of the mature trap. (AVI) S11 Movie. Integrated model colour-coded for cell area (rescale). Cellular-level integrated directional and areal conflict model coloured for cell area as shown in S10 Movie but continually rescaled to show shape change normalised for size. (AVI) S12 Movie. Areal conflict model with cell area. Cellular-level areal conflict model coloured for cell area, grey region shows approximate location of mouth ( Fig 9E, S5A Fig, S5B Fig). Scale bar is held constant to show increase in size from initial to resultant shape. Final shape is rotated at the end of the movie. Note enlarged cells in the bulge of the ventral midline region of the mature trap. (AVI) S13 Movie. Directional conflict model with cell area. Cellular-level directional conflict model coloured for cell area; grey region shows approximate location of mouth ( Fig 9F, S5C  Fig, S5D Fig). Scale bar is held constant to show increase in size from initial to resultant shape. heat shock (MP4) S23 Movie. Assigning polarity to quadrifid glands. Point coordinates at the quadrifid gland centre and at ends of each quadrifid gland arm were placed in VolViewer. Arrowheads were assigned oriented toward the greatest distance between arms with quadrifidScript software (DistArms) (Fig 12A-12D). (MP4) S24 Movie. Quadrifid gland polarity trap side view. OPT volume view of a U. gibba trap clipped to view quadrifid glands. Arrows flow from stalk to mouth through the side of the trap. Lines with no arrowheads were allocated when the difference in distance between arms was less than a threshold value of 2 μm (Fig 12E and 12F). OPT, Optical Projection Tomography. Lines with no arrowheads were allocated when the difference in distance between arms was less than a threshold value of 2 μm (Fig 12E and 12F).  Table. Goldengate construction. Breakdown of Goldengate parts used to create constructs used. Level 0 parts were synthesised and used to create level 0.5 and 1 constructs, and these were combined to make level 2 constructs, as described in [69]. (XLSX)