Organ geometry channels reproductive cell fate in the Arabidopsis ovule primordium
Abstract
In multicellular organisms, sexual reproduction requires the separation of the germline from the soma. In flowering plants, the female germline precursor differentiates as a single spore mother cell (SMC) as the ovule primordium forms. Here, we explored how organ growth contributes to SMC differentiation. We generated 92 annotated 3D images at cellular resolution in Arabidopsis. We identified the spatiotemporal pattern of cell division that acts in a domainspecific manner as the primordium forms. Tissue growth models uncovered plausible morphogenetic principles involving a spatially confined growth signal, differential mechanical properties, and cell growth anisotropy. Our analysis revealed that SMC characteristics first arise in more than one cell but SMC fate becomes progressively restricted to a single cell during organ growth. Altered primordium geometry coincided with a delay in the fate restriction process in katanin mutants. Altogether, our study suggests that tissue geometry channels reproductive cell fate in the Arabidopsis ovule primordium.
Introduction
A hallmark of sexual reproduction in multicellular organisms is the separation of the germline from the soma. In animals, primordial germ cells (PGCs) are setaside during embryogenesis from a mass of pluripotent cells. The number of germ cells depends on the balance between proliferation (selfrenewal) and differentiation, a process controlled by both intrinsic factors and signals from the surrounding somatic tissues. In flowering plants, the first cells representing the germline, the spore mother cells (SMCs), differentiate only late in development. SMCs arise multiple times, in each flower during the formation of the reproductive organs. In Arabidopsis, the female SMC differentiates in the nucellus of the ovule primordium, a digitshaped organ that emerges from the placental tissue of the gynoecium. The SMC is recognizable as a single, large, and elongated subepidermal cell, which is centrally positioned within the nucellus and displays a prominent nucleus and nucleolus (Bajon et al., 1999; Bowman, 1993; Schmidt et al., 2015; Schneitz et al., 1995).
Although SMC singleness may appear to be robust, more than one SMC candidate per primordium is occasionally seen, yet at different frequencies depending on the specific Arabidopsis accession (~5% in Landsberg erecta [Ler], 10% in Columbia [Col0], 27% in Monterrosso [Mr0]), (Grossniklaus and Schneitz, 1998; RodríguezLeal et al., 2015). Arabidopsis, maize, and rice mutants in which SMC singleness is compromised have unveiled the role of regulatory pathways involving intercellular signaling, small RNAs, as well as DNA and histone methylation (GarciaAguilar et al., 2010; Mendes et al., 2020; Nonomura et al., 2003; OlmedoMonfil et al., 2010; Schmidt et al., 2011; Sheridan et al., 1996; Sheridan et al., 1999; Su et al., 2020; Su et al., 2017; Zhao et al., 2008). As the SMC forms, cellcycle regulation contributes to the stabilization of its fate in a cellautonomous manner through cyclindependent kinase (CDK) inhibitors and RETINOBLASTOMARELATED1 (RBR1) (Cao et al., 2018; Zhao et al., 2017). SMC singleness thus appears to result from a twostep control: first, by restricting differentiation to one cell and second, by preventing selfrenewal before meiosis (reviewed in Lora et al., 2019; Pinto et al., 2019).
However, the precise mechanisms underlying the plasticity in the number of SMC candidates and SMC specification are still poorly understood. In principle, SMC singleness may be controlled by successive molecular cues. However, even in that scenario, such cues must be positional, at least to some extent, and thus involve a spatial component. Over the last decade, many different molecular cues defining spatial patterns in the ovule primordium were identified (Pinto et al., 2019; Su et al., 2020); however, their coordination is unknown. Since SMCs emerge at the primordium apex concomitant with its elongation, we hypothesize that geometric constraints during ovule morphogenesis influence SMC singleness and differentiation. Such a hypothesis could explain variation in the number of SMC candidates, ultimately culminating in a single SMC entering meiosis. Answering the questions of whether SMC formation follows a stereotypical or plastic developmental process and whether it is intrinsically linked to or independent of ovule primordium formation would unravel fundamental principles connecting cell fate establishment and organ growth.
Such an analysis requires a highresolution description of ovule geometry during development. Our current knowledge of ovule primordium growth in Arabidopsis is based on twodimensional (2D) micrographs from tissue sections or clearing. It is described in discrete developmental stages capturing classes of primordia by their global shape and SMC appearance until meiosis and by the presence of integument layers and ovule curvature later on (Grossniklaus and Schneitz, 1998). In addition, a 3D analysis of average cell volumes during primordium growth was recently provided (Lora et al., 2017), and extensive 3D analysis was carried on for late ovule stages (Vijayan et al., 2021). Yet, we lack a view of the patterning processes regulating early ovule primordium formation and how the dynamics of cell proliferation contributes to the cellular organization during primordium growth. Thus, we described and quantified the growth of the ovule primordium at cellular resolution in 3D. We combined 3D imaging, quantitative analysis of cell and tissue characteristics, reporter gene analyses, and 2D mechanical growth simulations. In addition, using the katanin mutant that affects anisotropic cell growth and division patterns (Luptovčiak et al., 2017a; Ovečka et al., 2020), we show that altered ovule morphology leads to ectopic SMC candidates. We also uncovered that differentiation of SMC candidates initiate earlier than previously thought, and provide evidence for a gradual process of cell fate restriction, channeling the specification of a single SMC prior to meiosis.
Results
Building a reference image dataset capturing ovule primordium development at cellular resolution
To generate a reference image dataset describing ovule primordium development in 3D and with cellular resolution, we imaged primordia at consecutive stages in intact carpels by confocal microscopy. Carpels were cleared and stained for cell boundaries using a modified PSPI staining (Truernit et al., 2008) and mounted using a procedure preserving their 3D integrity (Mendocilla Sato and Baroux, 2017; Figure 1A). We selected high signaltonoise ratio images and segmented them based on cell boundary signals using Imaris (Bitplane, Switzerland) as described previously (Mendocilla Sato and Baroux, 2017; Figure 1B). We manually curated 92 ovules representing seven consecutive developmental stages (7–21 ovules per stage, Figure 1B, Table 1, Figure 1—source data 1) and classified them according to an extended nomenclature (explained in Materials and methods). The temporal resolution of our analysis led us to subdivide early stages (stage 0I to stage 0III) covering primordium emergence prior to the straight digitshape of the organ set as stage 1I, where the SMC becomes distinguishable by its apparent larger size in longitudinal views (Grossniklaus and Schneitz, 1998; Figure 1B).
To evaluate the distinct contribution of domain, layer, and cellspecific growth dynamics, we labeled cells depending on their location in different regions of the ovule primordium: apical vs basal domains, L1,L2,L3 layers. In addition, we associated each cell with a cell type: ‘L1 apical’, ‘L1 basal’, ‘L2,L3 apical’, ‘L2,L3 basal’, ‘SMC’, ‘L1 dome’ (for the upmost apical L1 cells in contact with the SMC), ‘CC’ (for companion cells, elongated L2 cells adjacent to the SMC) (Figure 1B, Materials and methods).
To generate a quantitative description of ovule primordia with respect to cell number, size, and shape according to cell labels, layers, domains, ovule stage, and genotype, we developed an interactive, Rbased interface named OvuleViz. The interface imports cell descriptors exported from segmented image files and enables multiple plots from a userbased selection of (sub)datasets (Figure 1—figure supplement 1A–C, Materials and methods). This work generated a reference collection of annotated, 3D images capturing ovule primordium development at cellular resolution from emergence until the onset of meiosis. The collection of 92 segmented images, comprising a total of 7763 annotated cells and five morphological cell descriptors (volume, area, sphericity, prolate and oblate ellipticity), provides a unique resource for morphodynamic analyses of ovule primordium growth.
To identify correlations between growth patterns and differentiation, we first performed a principal component analysis (PCA) based on the aforementioned cell descriptors, per cell type and stage, considered together or separately (Figure 1—figure supplement 1D–E). In this global analysis, the SMC appears morphologically distinct at late stages (2I and 2II) but not at early stages. This prompted us to investigate in detail the contribution of different layers, domains, and cell types to ovule primordium growth and in relation to SMC differentiation.
Ovule primordium morphogenesis involves domainspecific cell proliferation and anisotropic cell shape patterns
The ovule primordium emerges from the placenta as a small domeshaped protrusion and grows into a digitshaped primordium with nearly cylindrical symmetry (stage 1I) before enlarging at the base (Figure 1B). Using our segmented images, we first quantified global changes in cell number, cell volume, and ovule primordium shape. Our analysis revealed two distinct phases of morphological events. Phase I (stages 0I to 0III) is characterized by a 4.5fold increase in total cell number together with a moderate increase in mean cell volume (10%, p=0.03). By contrast, Phase II (stages 1I to 2II) is characterized by a moderate increase in cell number (50%) and the global mean cell volume is relatively constant (Figure 2A, Figure 2—figure supplement 1A). To quantify the resulting changes in organ shape, we extrapolated a continuous surface mesh of the ovule outline and used it to compute its height and width at the base (Figure 2B–C, Figure 2—figure supplement 1B, Appendix 1). Anisotropic organ growth during Phase I was confirmed by a steady increase in height, while primordium width increased moderately (Figure 2C). This contrast in events between Phase I and II is illustrated by the foldchanges (FCs) in cell number and aspect ratio (Figure 2D), which range between 1.5 and 2.0 in Phase I but drop to 1.4 and 1.2 in Phase II. These observations confirmed that Phase I shows distinct growth dynamics compared to Phase II.
Next, to capture possible specific patterns of growth, we analyzed cell number, cell size, and cell shape using different viewpoints: one comparing the L1 and L2L3 layers and one contrasting the apical vs. basal domains. Counting cell number per viewpoint clearly showed a dominant contribution of the epidermis (L1) relative to the subepidermal layers and of the basal relative to the apical domain (Figure 2E and F). To verify these findings with a cellular marker, we analyzed the Mphasespecific promCYCB1.1::CYCB1.1dbGFP reporter (abbreviated CYCB1.1dbGFP) (UbedaTomás et al., 2009). We scored the number of GFPexpressing cells among 481 ovules and plotted relative mitotic frequencies per cell layer and domain for each ovule stage to generate a cellbased mitotic activity map (Figure 2G–H, Figure 2—source data 1, Materials and methods). In this approach, subepidermal (L2) cells beneath the dome were distinguished from underlying L3 cells to gain resolution in the L2 apical domain where the SMC differentiates. Consistent with our previous observation, in Phase I, a high proliferation activity was scored in L1 cells at the primordium apex (scoring 64% of all mitotic events). By contrast, the L2 apical domain remains relatively quiescent (only 3% of the mitotic events). During Phase II, the majority (60%) of mitotic events is found in the basal domain, consistent with the progressive population of the basal domain with more cells. It is of note that during this phase, few mitotic events are detected in L2 apical cells, with the exception of SMC neighbor cells that show frequent divisions at stage 1II. Thus, reporter analysis confirmed a biphasic, temporal pattern of cell division with changing regional contributions, suggesting the L1 dome and the basal domain as consecutive sites of proliferation, contributing to the morphological changes in Phase I and II, respectively.
Average cell size analysis, by contrast, did not reveal significant changes during primordium development with the notable exception of the SMC (Figure 2A, Figure 2—figure supplement 1C). The distinct size of the SMC candidate is already detected at stage 0III when compared to other L2,L3 cells (Figure 2I), or even earlier (stage 0I) when compared to all other cells (Figure 2I inset). Size differentiation of the SMC is not uniform among ovules, demonstrating plasticity in the process (Figure 2—figure supplement 1D). In addition, cells from the L2 and L3 layers are larger than L1 cells already at stage 0I (Figure 2J, Figure 2—figure supplement 1E), possibly due to a longer growth phase, consistent with the low division frequency observed previously.
We then investigated cell shape changes during primordium elongation, using ellipticity and sphericity indices computed following segmentation. The analysis did not reveal significant differences between domains or layers (Figure 2—figure supplement 1F). This could indicate either a highly variable cell shape or local, cellspecific differences. We thus more specifically analyzed the subepidermal domain where the SMC differentiates. Companion cells showed an increasing ellipticity (and decreasing sphericity), starting at stage 1I and culminating at stage 2I (Figure 2—figure supplement 1F). By contrast, the SMC only showed a moderate decrease in sphericity at late stages and no distinctive ellipticity at early stages (Figure 2—figure supplement 1F) when compared to other cells. To get more information on SMC shape, we compared the maximum, medium, and minimum anisotropy index. For this, we developed an extension for MorphoMechanX (Barbier de Reuille et al., 2015) (http://www.morphomechanx.org), (i) to perform a semiautomatic labeling of cell layers and cell types from a cellularized mesh obtained from segmentation data, and (ii) to compute in each 3D cell the principal axes of shape anisotropy and the corresponding indices (Figure 2—figure supplement 1G, Appendix 1, Figure 2—video 1). The averaged maximum anisotropy shape index of the SMCs was consistently above the medium (and hence above the minimum) anisotropy index (Figure 2K, Figure 2—figure supplement 1H). We measured the degree of alignment of the SMC major axis with the main growth axis of the ovule during early stages (stage 0II) and found a mean angle of 22° (±11°, n = 16) (Figure 2K). This confirmed that the SMC has a consistent anisotropic shape from early stages onwards, with a distinguishable major axis aligned with the primordium axis.
Taken together, these results suggest that anisotropic primordium growth is linked to a biphasic, domainspecific cell proliferation pattern, alternating between the L1 dome at Phase I and the basal domain at Phase II, combined with localized, anisotropic expansion in the L2 apical domain. In this process, SMC characteristics, such as distinct size, anisotropic shape, and orientation aligned with the growth axis of the primordium, emerge already in Phase I. The pronounced growth and elongation of the SMC in Phase II then occurs concomitant with primordium elongation. While primordium elongation is not explained by anisotropic cell growth alone but also depends on cell proliferation as shown above, the observation that cells are elliptic suggests a potential role for anisotropic cell growth. We explore this property in the next section through an in silico approach.
2D mechanical simulations relate ovule primordium growth to SMC shape emergence
Organ shape is determined by the rate and direction of cell growth, which is affected by signaling and the mechanical state and geometry of the tissue. This provides room for multiple regulatory feedback mechanisms and interactions between them (Bassel et al., 2014; Echevin et al., 2019). Mechanical constraints arise from the growth process in form of tensile and compressive forces that, in turn, influence cell and tissue growth (Echevin et al., 2019). To determine the role of ovule primordium growth on SMC differentiation, we sought to understand the contributions of local growth rate and anisotropy, and their relation to signaling and mechanical constraints (Coen et al., 2004; Kennaway et al., 2011).
For this, we developed two complementary 2D models of a longitudinal section of the ovule intersecting its main elongation axis: (i) a finite element method FEMbased mechanical model of the ovule represented as a continuous object (only outlining L1 vs. inner L2,L3 tissue) and (ii) a mass spring MSbased model of the ovule able to represent the mechanical status of each cell wall.
While MSbased models allow the investigation of the connection between organ and individual cell growth, the FEMbased models allow testing the role of material anisotropy. In addition, the two methods represent the cell wall in complementary ways. In FEM models, the cell wall is a continuous material throughout the tissue representation, while in MS models, the cell wall is modeled as a network of connected elastoplastic wires. The two approaches together allow us to determine the most plausible morphogenetic principles of ovule primordium growth, while addressing the contribution of specific parameters, such as material and cellbased properties in FEM and MS models, respectively.
Growth was implemented in both frameworks using two uncoupled but complementary modalities: (i) growth in which an abstract growth factor captures the cumulative effects of biochemical signals without considering their explicit mode of action (i.e. cell wall loosening, increasing turgor pressure), hereafter referred to as ‘signalbased growth’ (Boudon et al., 2015; Coen et al., 2004); and (ii) passive growth through relaxation of excess of strain inside the tissue, hereafter referred as ‘strainbased growth’ (Boudon et al., 2015; Bozorg et al., 2016).
The mechanical equilibrium is computed to ensure compatibility within the tissue that is locally growing at different rates and orientations. As a consequence, residual internal compressions and tensions arise and this, in turn, affects the mechanical behavior (Boudon et al., 2015; Rodriguez et al., 1994). A polarization field is used to set the direction of anisotropic growth (Coen et al., 2004). The simulation consists of cyclic iterations where the mechanical equilibrium is computed before each growth step is specified by signalbased growth or strainbased growth (Appendix 1). This strategy extends previous tissue growth models (Bassel et al., 2014; Boudon et al., 2015; Bozorg et al., 2016; Kuchen et al., 2012; Mosca et al., 2018). In the case of MS models, a cell will keep growing until it reaches a preassigned, userdefined target area, after which it will divide (shortest wall through the centroid rule, see Appendix 1, Figure 3—figure supplement 1D).
We designed a starting template consisting of an L1 layer distinct from the underlying L2,L3 tissue based on different growth and material properties (Table 2, Figure 3—figure supplement 1, Appendix 1). We first set the model components to produce a realistic, elongated, digitshaped primordium with a narrow dome and an L1 layer of stable thickness during development, fitting experimental observations (Figure 2—figure supplement 1C) (Reference Model, FEMModel 1, MSModel 2) (Figure 3A, C and E, Figure 3—figure supplement 1A–D). This model combines the following hypotheses: an initial, narrow domain of anisotropic, signalbased growth with a high concentration in inner layers, a broad domain competent for passive strainbased relaxation, and material anisotropy (for FEM models only). In addition, in the MSbased model, growth was prescribed to occur exclusively along the periclinal direction in the L1, according to our observations.
Then, using the versatility of the modeling framework to vary initial conditions, we tested the influence of the spatial distribution of the specified growth signal on primordium growth. As first variation, we let the growth signal diffuse broadly in the domain while maintaining the selected initial cells at a prescribed high intensity of the growth signal (FEMModel 4, MSModel 4). The emerging primordium is appreciably broader than the Reference Model (Figure 3—figure supplement 1B). We then explored the contribution of the growth signal in the L1 compared to the inner L2,L3 tissue to primordium growth. In FEMModel 5 and MSModel 5, the prescribed high growth signal is present only in the L1 but is free to diffuse to the L2,L3 layers. This produced a sharp primordium, narrower and taller than the primordium in the Reference Model, and a thicker L1 layer in the FEM model (Figure 3—figure supplement 1B, Figure 3—video 1). The L2 apical domain is narrower as compared to the Reference Models. When the growth signal is absent in the inner L2,L3 tissue (FEMModel 6, MSModel 6), signal growth in the L1 alone is not sufficient to enable primordium growth (Figure 3—figure supplement 1B).
To answer the complementary question whether a growth signal is required in the L1 for primordium growth when it is present in the L2,L3 layers, it was removed in FEMModel 7 and MSModel 7 (Table 2; Figure 3—figure supplement 1B). This scenario indeed enables growth of a digitshaped structure (as long as strainbased growth is permitted), yet the dome appears shallower than in the Reference Model in FEM simulations while it is comparable in MS simulations.
To conclude, both models where a high level of growth signal is selectively present in L1 or in inner L2,L3 layers can produce a digitshaped primordium. Yet, absence of a growth signal in the inner L2,L3 layers results in drastic shape alterations. This favors a scenario where inner L2,L3 layerdriven growth plays a fundamental role.
Next, through modulation of passive strainbased growth, we determined that a broad tissue domain uniformly competent for strain accommodation is necessary to resolve the high accumulation of stress, which limits primordium elongation (FEMModel 8, MSModel 8, Figure 3—figure supplement 1C). We also explored the contribution of material anisotropy for FEMbased simulations (FEMModel 2 and FEMModel 2a, Figure 3—figure supplement 1C): even in the case of isotropic material, it is possible to grow a digitshaped protrusion, yet with a wider dome and increased L1 thickness, contrasting experimental observations. To restore L1 thickness, it is sufficient to prescribe material anisotropy exclusively in the L1 (FEMModel 2a).
Finally, we asked whether growth anisotropy must be specified in the model or if organ geometry and mechanical constraints are sufficient to specify primordium shape. Removing the growth anisotropy component abolished the digit shape of the primordium and produced a hemispherical protrusion in both FEM and MSbased models (Figure 3B–F, Figure 3—video 2).
In summary, we identified parsimonious growth principles shaping the ovule primordium and suggesting different contributions of the epidermis and inner layers: an active tissue growth, mostly inner L2,L3 layerdriven and requiring a narrow, pitshaped domain of a growthsignal, complemented by passive tissue growth with a necessary response of the L1 to accommodate the accumulated strain. Furthermore, material anisotropy in the L1 is predicted to play a role in constraining L1 thickness as observed experimentally. The fact that two distinct modeling approaches converge on the same morphogenetic principles suggests that the proposed growth mechanisms are robust. Furthermore, when compared to the growth dynamics of real ovules from phase 0I to phase 1I, the cellbased MS model showed a good agreement (Figure 3—figure supplement 1D).
Next, we used cellbased MS simulations to explore the correlation between organ growth and SMC morphological differentiation. In the Reference Model (MSModel 2), an L2 apical cell with a trapezoidal shape, elongated along the main direction of ovule growth, emerged consistently during simulation (note that these cells still divide in the models as no special rule has been assigned to them) (Figure 3D). These are similar to SMC candidates at stage 0II in real primordia, in a 2D longitudinal, median section through the ovule (Figure 2K). The elongatedtrapezoidal shape of such a cell is not a prescribed feature of the model, but rather emerges from the combination of assigned anisotropic cell growth and geometrical constraints imposed by the surrounding, growing tissues.
Next, we explored the role of cell growth anisotropy in ovule primordium shape and SMC emergence. The MSModel 3 (Figure 3E–F) corresponds to a virtual mutant where cell growth is isotropic in inner layers (the L1 maintained anisotropic growth to preserve its thickness). This led to an ovule primordium with a wider and flatter dome, comparable to FEMModel 3. Despite the absence of a specified growth direction in inner L2,L3 tissue, due to geometrical constraints, the primordium still grows mostly vertically. We wanted to assess whether such geometrical constraints enable the formation of an elongated, trapezoidal SMC candidate also in the case of prescribed isotropic growth. As Figure 3G–H shows, the SMC candidate does not display the stereotypical shape and is not even elongated. Yet, mild cell anisotropy can be reached if the SMC candidate is allowed to grow twice more than in the Reference Model (MSModel 3a, Figure 3—figure supplement 1E).
Altogether, the different simulations suggest that the anisotropy and characteristic shape of the SMC could be an emerging property of ovule primordium growth connected to geometrical constraints, even in the absence of specified anisotropic growth. The complementary model, altering cell growth anisotropy specifically in the L1, further suggested that directional cell growth in the epidermis is necessary to accommodate inner L2,L3 layerdriven growth and permit primordium elongation (MSModel 3b, Figure 3—figure supplement 1F).
In all simulations, the candidate SMC eventually divided as the models miss a causative rule to differentially regulate cell division. When we prevented the SMC candidate to divide, its enlargement overrode primordium shape control in our simulations, creating an enlarged dome (Figure 3—figure supplement 1G). This suggests that a mechanism may limit SMC growth in real ovule primordia.
Altogether, 2D mechanical simulations in both continuous tissuebased (FEM) or cellbased (MS) approaches confirmed a role for localized cell growth compatible with experimental observations. The simulations also pointed to the importance of a differential role for the epidermis (accommodation) and the inner L2,L3 tissue (major growth component) as well as anisotropic cell growth as necessary for ovule primordium shape. Furthermore, the simulations suggested that SMC formation can be an emerging property of primordium geometry.
Katanin mutants show a distinct ovule primordium geometry
To experimentally test the prediction of the isotropic growth models, we analyzed ovule primordium growth and SMC fate establishment in katanin (kat) mutants with welldescribed and understood geometric defects: in absence of the microtubulesevering protein KATANIN, the selforganization of cortical microtubules in parallel arrays is hindered, thereby decreasing the cellulosedependent mechanical anisotropy of the cell wall, resulting in more isotropic growth (Bichet et al., 2001; Burk and Ye, 2002). Note that KAT is expressed ubiquitously and, thus, also in the ovule: the KATANIN protein could be detected in both epidermal and internal L2,L3 tissue layers using a GFP reporter (Figure 4—figure supplement 1).
Here, we specifically analyzed the shape and cellular organization of ovule primordia in kat mutants using the botero (bot17) (Ws background), lue1, and mad5 alleles (Col background) (Bichet et al., 2001; Bouquin et al., 2003; Brodersen et al., 2008). We generated and analyzed a new dataset of 59 annotated 3D digital kat and corresponding wildtype ovule primordia at stages 0III, 1I, and 1II (Figure 4—source data 1). kat mutant primordia clearly showed an increased size and a more isotropic shape (Figure 4A), being 1.5 times bigger in volume than wildtype primordia (p=0.007, stage 0III, Figure 4B) with a smaller aspect ratio (p<0.01 stages 1I, 1II, Figure 4C). Because the widthtoheight ratio does not inform on the shape at the flanks, we derived an equation to estimate the ‘plumpiness’ of the primordia: primordia with rounder flanks will have a higher bounding box occupancy, that is, the volume fraction of a fitting, 3D parallelepiped (bounding box) effectively occupied by the primordium (Figure 4D, left), than straight digitshaped ovules of the same aspect ratio. Mutant primordia are clearly distinct from wildtype primordia in their relationship between aspect ratio and bounding box occupancy at stage 1I (Figure 4D, Figure 4—figure supplement 2A), when the primordium normally starts elongating along the major growth axis. These measurements confirmed a marked attenuation of anisotropic growth in kat ovule primordia as was observed in roots, shoot organs, or seeds (Bichet et al., 2001; Hervieux et al., 2016; Luptovčiak et al., 2017b; Ren et al., 2017; Uyttewaal et al., 2012; Wightman et al., 2013; Zhang et al., 2013). Yet, at a global level, the mean cell number, cell volume, and sphericity are not significantly different between kat and wildtype primordia (Figure 4—figure supplement 2B–C). Thus, these global approaches did not resolve the origin of primordia misshaping.
To refine the analysis, we contrasted different layers as done for wildtype primordia and found local alterations that provide an explanation for the formation of broader and more isotropic primordia in kat mutants. Indeed, at stage 0III, kat ovules display a broader epidermis domain composed of more and larger cells than wildtype ovules (p=0.03 and p<0.001, respectively) (Figure 4E). Mitosis frequency analysis indicated a shift in cell division from the apex toward the basis (2.4 times less mitoses in the L1 dome domain and 2.7 times more in the L1 basal domain, compared to the wild type) (Figure 4F, Figure 4—figure supplement 2E), consistent with the increased cell number observed in L1. Yet increased cell division is not a general characteristic of kat primordia since, overall, the apical and basal domains are not massively overpopulated. By contrast, kat cells are generally larger and slightly more spherical in all domains (Figure 4E, Figure 4—figure supplement 2C–D). When specifically looking at the L2 apical domain, where the SMC differentiates, we noticed an increased relative frequency of CYCB1.1dbGFP expression in kat SMC neighbors at stage 0II, whereas it decreases at stage 1II (Figure 4F, Figure 4—figure supplement 2E). This is in stark contrast with SMC neighbors in wildtype primordia displaying first a relative mitotic quiescence at stage 0II, and then enhanced mitotic activity at stage 1II. Yet, in the SMC, no mitotic activity increase was observed in kat as compared to wildtype primordia (Figure 4F, Figure 4—figure supplement 2E).
In conclusion, the absence of KATmediated cell growth anisotropy is associated with spatiotemporal shifts in cell divisions, leading to an altered primordium geometry including a flatter dome and enlarged basis. Interestingly, these changes coincided with the occurrence of additional large cells in the L2 apical domain of kat primordia (SMC direct neighbor cells, Figure 4G, Materials and methods) similar in size to the SMC. This raises the question of the identity of these ectopic, enlarged neighbor cells.
Altered ovule primordium geometry in katanin mutants induces ectopic SMC fate
Following the observation of misshaped kat primordia that contained larger L2 apical cells, we asked whether this had an impact on cell identity and SMC establishment. As in the wild type, a clear SMC is identifiable in kat primordia, yet it appears slightly bigger and more spherical with some variability over stages (Figure 5A, Figure 5—figure supplement 1). Interestingly, and consistent with the increased mitotic frequency in SMC neighbors observed at earlier stages, at stage 1II, we scored additional SMC neighbors in kat primordia. On average, these were larger (14%, p=0.006) and slightly more isotropic in shape (~11% less ellipsoid p=0.03, ~2% more spherical, p<0.01) than in the wild type (Figure 5B). At stage 1II, 34% of the kat primordia (n = 109) showed more than one enlarged, subepidermal cell; decreasing to 14% at stage 2II (n = 104), in stark contrast with wildtype primordia showing a majority (86% to 96%) of primordia with an unambiguous, single SMC (Figure 5C). The kat phenotype is thus reminiscent of mutants affecting SMC singleness (Pinto et al., 2019) and we hypothesized that these enlarged cells are ectopic SMC candidates.
To verify this hypothesis, we introgressed several markers of SMC identity in a kat mutant background. The first marker is a GFPtagged linker histone variant (H1.1GFP) that marks the somatictoreproductive fate transition by its eviction in the SMC at stage 1I (She et al., 2013). H1.1GFP eviction occurred in more than one cell in 29% of kat primordia (n = 43, Figure 5D). The second marker reports expression of the KNUCKLES transcription factor (KNUYFP) in the SMC (Tucker et al., 2012). Detectable as early as stage 1I in the wild type, it was ectopically expressed in 21% of kat primordia at stage 1II (n = 43) (Figure 5E, Figure 5—figure supplement 2A). Third, to test whether the ectopic SMC candidates have a meiotic competence, we analyzed the AtDMC1GUS reporter (Agashe et al., 2002; Klimyuk and Jones, 1997) and, indeed, scored ectopic expression in 57% of kat primordia (n = 56) (Figure 5F). However, using aniline blue staining of callose deposition, which in wildtype SMCs stains the cell wall immediately before meiosis and marks the cells walls of tetrads after meiosis, we never detected ectopic cells accumulating callose or ectopic tetrads in kat ovules (n = 119 and 269 ovules in bot17 and lue1, respectively) (Figure 5—figure supplement 2B). We also noted a significant frequency of kat ovules lacking callose in SMCs or tetrads in both kat alleles, suggesting altered cell wall composition in the SMC, consistent with previous reports on kat somatic tissues (Burk and Ye, 2002) We next analyzed the pWOX2CenH3GFP marker labeling centromeres starting at the functional megaspore stage (FG1 stage) (De Storme et al., 2016) and found 21% kat primordia (n = 29) with two labeled cells, as compared to 4% (n = 26) in the wild type (Figure 5—figure supplement 2C). These ectopic spores are haploid since they display five centromeres. They may correspond to ectopic surviving spores given their alignment with the basalmost functional megaspore (Figure 5—figure supplement 2C). Using an additional gametophytic reporter, pAKV:H2BYFP (Schmidt et al., 2011), we confirmed a significant number of ectopic spores at stage FG1 (31% in kat vs. 11% in wildtype ovules) showing a residual signal (Figure 5—figure supplement 2D). However, at later stages (FG2FG7), we did not find evidence for ectopic embryo sacs (Figure 5—figure supplement 2D). At the FG7 stage, notable alterations in ovule morphology and the number of nuclei in the embryo sac were visible in kat as previously reported (Luptovčiak et al., 2017b).
Taken together, these results show that ectopic, abnormally enlarged SMC neighbors in kat primordia show at least some characteristics of SMC identity, yet they do not complete meiosis and are likely reincorporated into the soma. After meiosis, ectopic surviving spores are observed but do not complete gametogenesis. Hence, reproductive fate establishment is altered in kat ovules but the defects do not persist beyond meiosis suggesting a regulative process.
SMC singleness is progressively resolved during primordium growth
Based on our analysis of kat primordia, it appeared that the frequency of ectopic SMC candidates was high at early stages but decreased over time (Figure 5C). This is reminiscent of the phenotypic plasticity in SMC differentiation observed, to a lesser degree, in different Arabidopsis accessions (RodríguezLeal et al., 2015). To characterize this plasticity during development, we analyzed 1276 wildtype primordia from stage 0II to 2II and scored the number of primordia with one or two enlarged, centrally positioned, subepidermal cells (Classes A and B, respectively, Figure 6A–B). In the wild type, the majority of ovules showed a single candidate SMC (Class A) at stage 0II but 27% of primordia (n = 289) had two SMC candidates (Class B), this frequency decreasing to 3% at stage 2II (n = 103). This finding is consistent with ~ 5% primordia in wildtype at stage 1II showing H1.1GFP eviction (n = 43) and KNUYFP expression (n = 88) in more than one cell, respectively (Figure 5D–E). Thus, instead of being immediate, SMC singleness can arise from a progressive restriction of fate among several SMC candidates during development. In wildtype primordia, SMC singleness is largely resolved at stage 1I (Figure 6B).
Next, we quantified Class A and B primordia in kat mutants by scoring 2587 ovules of three different mutant alleles (Figure 6C, Figure 6—figure supplement 1A). Clearly, plasticity is higher with 42% (n = 202) class B primordia at stage 0II in kat compared to 26% (n = 289) in the wild type. In addition, the resolution process is delayed in kat mutants as up to 33% Class B primordia persist at stage 1II (n = 113) compared to 12% in the wild type (n = 281). Consistently, the two SMC markers used previously, clearly identify ectopic SMC candidates in kat primordia at stage 1I (Figure 6—figure supplement 1B,C), and more significantly at stage 1II (28.5% ectopic eviction of H1.1GFP and 21% ectopic expression of KNUYFP, Figure 5D–E). In addition, at stage 2II, when 17% of cleared kat primordia (n = 110) showed ectopic SMCs, H1.1GFP eviction and KNUYFP expression was only found in 8% (n = 25) and 7% (n = 74) of the primordia, respectively (Figure 6—figure supplement 1B,C). Therefore, molecular events associated with SMC fate are partially uncoupled from cell growth during the resolution of SMC singleness in kat.
Another outcome of this study is the early emergence of SMC candidates, based on cell size mostly and confirming our former 3D cellbased analysis that showed increased cell volume of the L2,L3 apical cells already at stages 0I/0II (Figure 2J). To corroborate this finding using a different criterion, we measured the nuclear area, which is a distinctive feature of SMCs (RodríguezLeal et al., 2015; She et al., 2013). We compared the nuclear area of SMC candidates to that of surrounding L1 and L2 cells, in both Class A and B primordia, at stages 0II and 0III. Wildtype Class A and B primordia showed an enlarged nucleus in the candidate SMCs from stage 0II onwards, and this correlation was also true in kat primordia (Figure 6D, Figure 6—figure supplement 1D,E).
However, it remains difficult to resolve the precise timing of SMC establishment at these early stages due to the lack of molecular markers. Yet, we reasoned that we may be able to distinguish the SMC candidates from their neighbors by their cell cycle pattern, where cells entering meiosis may engage in a specific Sphase compared to regularly cycling mitotic cells. To this aim, we used a GFPtagged PCNA variant marking the replication machinery, pPCNA1::PCNA1:sGFP (PCNAGFP) (Yokoyama et al., 2016). When engaged at active replication forks, PCNAGFP shows nuclear speckles characteristic of Sphase (Strzalka and Ziemienowicz, 2011). During G1/G2, it remains in the nucleoplasm and in Mphase it is undetectable. We specifically quantified the distribution patterns of PCNAGFP in cells from the L2 apical domain in wildtype and kat primordia, separately for Class A and B ovules (Figure 6E–H).
In Phase I wildtype primordia, PCNAGFP was always detectable, in both classes, indicating that L2 apical cells are rarely in Mphase, consistent with the seldom detection of mitoses using CYCB1.1dbGFP. We observed an Sphase pattern consistently in one (Class A) or more (Class B), centrally positioned L2 apical cells, presumably corresponding to SMC candidates (Figure 6E). This pattern was captured in a large proportion of primordia at stage 0II and 0III: ~24% (n = 112) and ~40% (n = 79) for Class A, respectively; 30% (n = 41) and 28% (n = 60) for Class B, respectively (Figure 6F). Such high frequencies could be generated either by a slow Sphase in SMC candidates only (the persistence of the marker increasing the probability to score it repeatedly in our sample size), or by a regular (short) mitotic Sphase in a high number of SMC candidates. The low detection frequency of the mitotic marker CYCB1.1dbGFP in SMC candidates is inconsistent with the latter possibility. Thus, the likeliest interpretation is that candidate SMCs enter a slow Sphase from early stages onwards, probably meiotic, although this cannot be assessed with this marker.
In kat primordia at Phase I, Class A ovules showed a wildtype pattern in the SMC at stage 0II, but a reduction at stage 0III, suggesting alterations in Sphase entry; and Class B ovules, by contrast, displayed high frequency of Sphase patterns, indicating either a longer Sphase or multiple dividing cells (Figure 6—figure supplement 1G).
In Phase II, SMC candidates showed no Sphase pattern in a large majority of both Class A and B primordia (97% on average, over entire Phase II) (Figure 6G–H, Figure 6—figure supplement 1G–H) in agreement with the presence of newly replicated DNA at stage 1I (She et al., 2013). In contrast to Phase I, however, Sphase is now detected in SMC neighbors (~21% at stage 1II) (Figure 6G–H), consistent with the divisions observed in these cells (Figure 2D). Strikingly, in Phase II, kat primordia display a higher frequency of Sphase patterns in SMC candidates (Figure 6H, Figure 6—figure supplement 1F,H), suggesting that Sphase duration or entry is delayed compared to the wild type. SMC neighbors, by contrast, show a lower frequency of Sphase patterns (Class A), consistent with reduced division in these cells (Figure 4F).
Collectively, our data indicate a cellular heterogeneity in terms of size, nuclear size, and Sphase patterns, of the L2 apical domain as compared to L1, which leads to the emergence of one or several SMC candidates as early as stage 0II. The gradual decrease in the number of primordia with ambiguous SMC candidates demonstrates a developmentally regulated resolution of SMC fate to a single cell. This process is associated with a specific cell cycle progression, cellular elongation, and a robust expression of SMC fate markers. In kat primordia displaying alterations in geometry, SMC singleness is largely compromised: plasticity in SMC emergence is increased and fate resolution to a single SMC is delayed.
Discussion
Organogenesis involves coordinated cell division and cell expansion, complex growth processes orchestrated by biochemical and mechanical cues (Echevin et al., 2019). How cell differentiation is coordinated in space and time during organ growth, and whether these processes are interrelated, are central aspects for the elucidation of patterning principles (Whitewoods and Coen, 2017). The female germline is initiated with SMC differentiation in the ovule primordium. The SMC emerges as a large, elongated, subepidermal cell that is centrally located at the apex of the primordium, a digitshaped organ emerging from the placenta. To study how SMC fate relates to ovule organogenesis, we generated a reference collection of images capturing ovule primordium development at cellular resolution in 3D and determined cell division frequencies in the different domains. We observed a biphasic pattern of cell divisions alternating in the epidermis and inner layers, as well as the apical and basal domains, in Phase I (stages 0I to 0III) and Phase II (stages 1I to 2II).
However, this approach did not allow us to resolve the driving morphogenetic factors. For this reason, we developed continuous tissuebased and cellbased 2D simulations of primordium growth. The different simulations revealed key growth principles shaping the ovule primordium and uncovered differential roles for the epidermis and inner layers. Notably, an inner tissuedriven growth model, where the L1 also contributes to the expansion of the primordium, best described ovule primordium growth. This is reminiscent of a model describing leaf primordium emergence (Peaucelle et al., 2011). In addition, bestfit models produced by both cell and tissuebased simulations predicted a growthpromoting signal in a confined domain along a vertical stripe at primordium emergence. Candidate growth signals are phytohormones, peptides, and small RNAs known to affect ovule primordium growth (Kawamoto et al., 2020; Pinto et al., 2019; Su et al., 2020). The domains of auxin response restricted in the L1 dome and of cytokinin signaling localized in a region basal to the SMC in Phase II primordia (Bencivenga et al., 2012) also suggest a confined growth signal. Whether the signaling domains are established already in Phase I and play a causative role in primordium patterning remains to be determined. The epidermis, by contrast, is predicted to play a key role in accommodating the constraints generated by innertissue growth. In this layer, strainbased growth and anisotropic material properties possibly resolve mechanical conflicts arising between tissue layers that grow at different rates (Hervieux et al., 2017). In line with this hypothesis, we observed frequent divisions in L1 apical cells in vivo that support the expansion of the epidermis while inner tissues develop.
Interestingly, while our models did initially not contain an a priori rule to produce the typical, elongated shape of the SMC, it emerged consistently as a trapezoidshaped L2 cell in the cellbased reference model. This shape likely emerges from the combination of assigned anisotropic cell growth and geometrical constraints imposed by the surrounding growing tissues. Explicitly blocking SMC division during the simulation, did not only enable its expansion as expected, but also pushed surrounding cells and strongly deformed ovule morphology. Thus, ovule growth homeostasis in vivo likely requires a mechanism to accommodate the differential growth of the SMC. A gradual reduction of SMC turgor pressure is a plausible scenario that would limit SMC size and prevent overriding the constraints provided by surrounding cells, similar to what was suggested for the shoot apical meristem (Long et al., 2020). In turn, a gradient of pressure in a field of cells could provide positional information through the directional movement of water and other molecules, thereby linking organ growth homeostasis, cell growth, and cell fate (Beauzamy et al., 2014; Long et al., 2020). Such a mechanism could also participate in determining the domain of the growth signal predicted by our simulations.
Another prediction of our growth models is the key role of anisotropic cell growth in controlling the geometry of the ovule primordium. Primordia of the kat mutant, deficient in the microtubule severing enzyme KATANIN (Luptovčiak et al., 2017a), resemble virtual mutant primordia generated by models where cell growth is isotropic in the inner layers. kat primordia have a flatter dome and large basis associated with global alterations in the cell proliferation pattern, cell size, and cell shape. Interestingly, kat primordia develop ectopic SMC candidates as early as stage 0II. The most parsimonious hypothesis is that the altered geometry in kat primordia expands the domain of cells competent to form SMCs. Yet, we cannot exclude a direct effect of the kat mutation on L2 apical cells, disconnected from organ geometry, which would induce de novo SMC fate. However, this scenario is unlikely. First, KATANIN is present throughout ovule primordium cells. Second, we would expect increased cell growth and slower mitoses (Luptovčiak et al., 2017a), resulting in a reduced division frequency of L2 apical cells, which is not the case. Instead, we measured increased cell divisions in L2 SMC neighbors at stage 0II. Also, the delayed divisions of SMC neighbors in kat at stage 1II cannot explain the formation of ectopic SMC candidates at stage 0II.
Therefore, the most parsimonious explanation is that the emergence of several SMC candidates in kat primordia is an effect of ovule primordium geometry. In this scenario, isotropic cell growth and altered mechanical constraints in the tissue acting from primordia emergence onwards, lead to divisions and patterning alterations that expand the domain of cells competent for SMC fate. This working model paves the way to explore the role of epidermal geometry in controlling regulators of the cell cycle and SMC fate in L2 apical cells. This is reminiscent of the epidermis in the shoot apical meristem, affecting the dome shape and acting on stem cell regulators in underlying layers (Gruel et al., 2016; SavaldiGoldstein et al., 2007). The predicted role of the primordium epidermis to accommodate – and perhaps feedback on – underlying growth constraints is particularly interesting considering the role of mechanical cues on gene regulation (Fal et al., 2017; Landrein et al., 2015). The epidermis of the ovule primordium is also a known source of signaling cues (Kawamoto et al., 2020; Pinto et al., 2019; Su et al., 2020). In addition, phytohormones act themselves on KATANINmediated oriented cell growth and cell division (Luptovčiak et al., 2017a). In this context, it is possible that misshaping of kat primordia arises from a disrupted feedback altering the distribution pattern of the hypothesized growth signals.
Furthermore, our analyses unveiled new characteristics of the SMC establishment process. We found that SMC candidates emerge from within a mitotically quiescent L2 apical domain, consistent with the finding that the H3.1 histone variant HTR13 is evicted, marking cell cycle exit (HernandezLagana and Autran, 2020). In addition, SMC candidates have a markedly long Sphase compared to surrounding cells. These observations are reminiscent of the animal germline where mitotic quiescence is a prerequisite for meiosis (Kimble, 2011; Reik and Surani, 2015). Also, early SMC candidates already display a typically large cell and nucleus size and an elongated shape aligned with the primordium growth axis. Collectively, we found that SMC characteristics are established much earlier than previously thought, that is, soon after primordium emergence (stages 0II/0 III). Moreover, these characteristics frequently arise in more than one SMC candidate at Phase I, typically resolving into a single SMC at the onset of Phase II. This clearly documents a developmental sequence of plasticity at SMC fate emergence and progressive resolution of SMC fate to a single cell. This is reminiscent of developmental canalization, which refers to the capacity of an organism to follow a given developmental trajectory in spite of disturbances (Hallgrimsson et al., 2019; Scharloo, 1991; Waddington, 1942). Cell fate canalization is well studied in animal systems where it is modulated by intercellular signalbased feedback mechanisms (Heitzler and Simpson, 1991), epigenetic regulation (Pujadas and Feinberg, 2012), and organ geometry (Huang et al., 2020; Huang and Russell, 1992; Pujadas and Feinberg, 2012; Royer et al., 2020). In plants, canalization is better known in the context of organogenesis during phyllotaxis and developmental robustness (Godin et al., 2020; Lempe et al., 2013). Our study expands the examples of canalization to the level of a cellular domain in the Arabidopsis ovule primordium. While L2 apical cells initially share the competence to form SMC candidates, leading to plasticity at SMC emergence, the progressive restriction of cell fate possibilities in the primordium apex ultimately leads to the specification of only one SMC committed to meiosis. Our results are in line with a formerly proposed canalization process operating during SMC establishment (Grossniklaus and Schneitz, 1998; RodríguezLeal et al., 2015). Despite the fact that several mutations (reviewed in Pinto et al., 2019; Mendes et al., 2020), including kat (this study), alter SMC singleness in Arabidopsis, canalization remains a robust process securing the formation of a single embryo sac despite such genetic perturbations. How this developmental mechanism buffers phenotypic interindividual variations and whether it is evolutionary constrained remains to be determined.
In this study, we quantified plasticity among ovule primordia and progressive fate resolution during primordium growth in wildtype reference accessions. We characterized specific cellular events associated with these processes, notably differential cell growth and cell division. SMC fate emergence is characterized by mitotic quiescence and cellular growth in one or more L2 apical cells. SMC singleness resolution is associated with reentry into a somatic cell cycle (this study), and reincorporation of a replicative histone H3.1 (HernandezLagana and Autran, 2020) in candidates neighboring the SMC. How known epigenetic and signaling factors interplay to secure SMC singleness remains to be determined. Similar to mouse and Drosophila, where tissue mechanics and organ geometry were shown to contribute to cell fate canalization (Chan et al., 2017; Huang et al., 2020; Royer et al., 2020), we propose that ovule primordium geometry contributes to channel SMC fate in the apex and the resolution into a single SMC. In this conceptual framework, kat increases plasticity and delays the resolution process toward SMC singleness (working model Figure 6I).
Altogether, our work proposes a conceptual framework linking organ geometry, cell shape, and cell fate acquisition in the ovule primordium, which is potentially of broader relevance in plant patterning. In addition, the image resource published in this study is complementary to others capturing ovule development at later stages (Lora et al., 2017; Vijayan et al., 2021). It also populates a growing number of 3Dsegmented images of plant tissues and organs (Wolny et al., 2020), which collectively build the fundament of a developmental atlas integrating morphogenesis with gene expression (Hartmann et al., 2020).
Materials and methods
Plant growth and plant material
Request a detailed protocolArabidopsis thaliana plants were grown under longday conditions (16 hr light) at 20–23°C in a plant growth room. Columbia (Col0) and Wassileskija (Ws4) accessions were used as wildtype controls depending on the mutant background used in the experiment. Three katanin alleles were used: bot17 (Bichet et al., 2001) in Ws4 accession, lue1 (Bouquin et al., 2003) and mad5 (Brodersen et al., 2008) both in the Columbia (Col0) accession. Homozygous mutant individuals for all the katanin alleles were identified on the basis of their recessive vegetative phenotype. The following published markers were used: pCYCB1.1:dbGFP (UbedaTomás et al., 2009), pKNU:nls:YFP (Tucker et al., 2012), pH1.1:H1.1:GFP (She et al., 2013), AtPCNA1:sGFP (Yokoyama et al., 2016), pAtDMC1:GUS (Agashe et al., 2002; Klimyuk and Jones, 1997), pWOX2:CenH3:GFP (De Storme et al., 2016), pAKV:H2B:GFP (Schmidt et al., 2011), and crossed to kat mutants, and to Ws4 ecotype for bot17 allele comparisons. For KATANIN localization, the published reporter line GFP:KTN1 in a ktn12 background (Lindeboom et al., 2013) was used.
3D imaging and image processing (segmentation and labeling)
Request a detailed protocolEntire carpels were stained using the pseudoSchiff propidium iodide (PSPI) cell wall staining procedure providing excellent optical transparency for 3D imaging in depth in wholemount. We described previously the manipulation, staining, mounting of the flower carpels and imaging procedures (Mendocilla Sato and Baroux, 2017). Cellboundary based image segmentation was done using ImarisCell (Bitplane) as described in details previously (Mendocilla Sato and Baroux, 2017). Each ovule was manually labeled in Imaris using customized Cell Labels for the different cell types and domains colored as shown in Figure 1. We defined the labels as follows:
SMC (Spore Mother Cell): most apical central enlarged L2 cell. At stage 0I, as enlargement is not always detected visually, the most apical L2 cell was then selected as candidate SMC (cSMC).
L1: epidermal cells
L1 dome: most apical cells in contact with SMC
L2,L3: cells below the epidermis. L2 and L3 were not distinguished originally.
Apical domain: group of cells at the apex of the primordium and encompassing the SMC and direct neighbor cells.
Basal domain: group of cells at the basis of the primordium below the apical domain and until, but not including cells of the placental surface. At stages 0I and 0II only an apical domain is defined. A basal domain appears only starting stage 0III.
CC (Companion Cell): L2 cells in apical domain in contact with the SMC with an elongated shape (as judged in ovule longitudinal median section using the ‘clipping plane’ IMARIS tool).
SMC contact: cells in contact with the SMC.
Semiautomated segmentation requiring user input and manual labeling can be error prone. To reduce the error rate, the 92 images were segmented and labeled by one author, but verified and curated by two others. All Imaris files used in the study are available at the DRYAD repository: https://doi.org/10.5061/dryad.02v6wwq2c.
Ovule stage classification
Request a detailed protocolOvule development is described according to a wellaccepted nomenclature (Schneitz et al., 1995). The first stage, initially defined as stage 1I, indistinctly grouped primordia from emergence until digit shape. To enable describing early morphogenetic processes, however, we propose to (i) restrict stage 1I to the final digit shape stage and (ii) subdivide preceding stages as stages 0I, 0II, and 0III, as shown in Figure 1B. These developmental stages are classified according to the approximate number of cell layers protruding above the placenta and overall shape of the ovule as described in Table 1 and Figure 1—source data 1.
Quantification of cell number, size and shape and interactive plotting using OvuleViz
Request a detailed protocolSeveral cell descriptors were retrieved using the ImarisCell’ Statistics function and exported as. csv files: cell area, cell volume, cell sphericity, cell ellipticity (oblate and prolate). We developed an interactive Rbased data plotting interface, OvuleViz, reading the Imarisderived data within the exported files ordered by genotype then stages (Figure 1—figure supplement 1). OvuleViz is freely available at https://github.com/barouxlab/OvuleViz (Pires, 2021), and is based on a shiny interface for R. OvuleViz allows plotting selectively one or several of the cell descriptors for chosen stages and genotypes, along different visualization (scatter plots, box plots, histograms). In addition, OvuleViz retrieves the cell number from the number of objects in the. csv file.
3D quantification of ovule volume and shape using IMARIS
Request a detailed protocolOvule volume and shape were quantified on 3D segmentations using IMARIS software, to compare katanin and wildtype genotypes. For each ovule, all segmented labeled cells were duplicated and fused as a single cell object. The ImarisCell’ Statistics function was used to retrieve ‘cell volume’ and ‘bounding box OO’ (objectoriented 3D bounding rectangle, exporting the minimum (A), mid (B) and maximum (C) lengths of the object bounding rectangle). Width to Height ratio was calculated by dividing the maximum (C) length by the mean of mid (B) + minimum (A) lengths. Bounding box occupancy was calculated as the ratio of ovule volume by the bounding box volume.
Modeling and image analysis with MorphoMechanX
Request a detailed protocolThe details of tissue growth models and MorphoMechanXbased processing are described in Appendix 1.
2D cytological analysis of cleared ovule primordia and quantifications
Request a detailed protocolFor cytological examination of cleared ovule primordia, flower buds from wildtype and mutant plants were harvested and fixed in formalinacetic acidalcohol solution (40% formaldehyde, glacial acetic acid, 50% ethanol; in a 5:5:90 vol ratio) for at least 24 hr at room temperature. After fixation, samples were washed two times with 100% ethanol and stored in 70% ethanol. Gynoecia of 0.2–0.6 mm in length were removed from the flowers with fine needles (1 mm insulin syringes), cleared in Herr’s solution (phenol: chloral hydrate: 85% lactic acid: xylene: clove oil in 1:1:1:0,5:1 proportions), and observed by differential interference contrast microscopy using a Zeiss Axioimager Z2 microscope and 40X or 63X oil immersion lenses. Picture acquisition was done with a sCMOS camera (Hamamatsu ORCA Flash V2). Nuclei area measurements were carried out with ImageJ software, using the manual contour tool ‘Oval’.
Fluorescence microscopy and quantifications
Request a detailed protocolEpifluorescence imaging of callose using aniline blue staining was performed as described (Cao et al., 2018), and observed on a Zeiss Axioimager Z2 microscope with CFP emission filter, DIC, and a 63X oil immersion lens. Image acquisition was done with a sCMOS camera (Hamamatsu ORCA Flash V2).
Imaging of ovule primordia stained in wholemount for cell boundary was done as described (Mendocilla Sato and Baroux, 2017) using a laser scanning confocal microscope Leica LCS SP8 equipped with a 63X glycerol immersion objective and HyD detectors.
Imaging of the GFP and YFP markers was performed using a laser scanning confocal microscope Leica LCS SP8 equipped with a 63X oil immersion objective, or a Leica SP5 equipped with a 63X glycerol immersion objective and HyD detectors. Samples were mounted in 5% glycerol with the cell wall dye Renaissance 2200 (SR2200) diluted 1/2000 or as described (Musielak et al., 2015). The following wavelengths were used for fluorescence excitation (exc) and detection of emission (em): Renaissance: exc = 405 nm, em = 415–476 nm; GFP: exc = 488 m, em = 493–550 nm; YFP: exc = 514 nm, em = 590–620 nm. For KATANIN localization, GFP:KTN1 ktn12 carpels were mounted in gelrite 0.2% supplemented with Renaissance (SR2200) diluted 1/1000 in water, and imaged immediately using a Zeiss LSM880 laser scanning microscope equipped with a 40X long distance water immersion objective and Airyscan detector, in the SR (superresolution) acquisition mode. Fluorescence was collected using the following filters settings: 405 nm and 488/561/633 nm primary beam splitters for all channels, then BP420−480 + BP495550 secondary filter for 488 nm excitation channel (GFP) and BP420−480 + LP605 for 405 nm excitation channel (Renaissance). Images were processed using the builtin Airyscan processing tool.
All images were processed using ImageJ (Rueden et al., 2017), OMERO (Besson et al., 2019), or Imaris (Bitplane, Switzerland) for contrast/intensity adjustments and maximum intensity projections where relevant.
The mitotic activity in both wildtype and katanin ovules was quantified by scoring the cells expressing promCYCB1::dbCYCB1GFP (M phase reporter) on 3D stacks, in each ovule domain at each developmental stage. Only ovules showing at least one cell expressing promCYCB1::dbCYCB1GFP were included in the analysis. At a given stage, the frequency of mitoses per domain is the ratio between the number of GFPpositive cells within a domain in all observed ovules, to the total number of GFPpositive cells found in all domains in that population of ovules.
To quantify ectopic expression of H1.1GFP, KNUnlsYFP, pWOX2:CenH3GFP and pAKV:H2BGFP markers, the percentage of ovules presenting expression (or eviction in the case of H1.1GFP) of the marker in more than one cell (or in more than one embryo sac for pAKV:H2BGFP) was scored on 3D stacks. To quantify PCNAGFP patterns, using 3D stacks, cells of the L2 apical domain were classified according to the ‘speckled’ or ‘nucleoplasmic’ patterns, or absence of the marker.
Histochemical detection of uidA reporter gene product (GUS staining)
Request a detailed protocolGynoecia of 0.4–0.6 mm in length were removed from the flowers with fine needles and placed in staining solution, using high stringency conditions for ferro and ferricyanide concentrations to limit GUS product diffusion (0.1% Triton X100, 10 mM EDTA, 5 mM ferrocyanide, 5 mM ferricyanide and 20 mg/ml 5 bromo4chloro3indolylbetadglucuronic acid cyclohexylammonium salt (Xgluc, Biosynth AG, Staad, CH) in 50 mM phosphate buffer), for 96 hr at 37°C. After staining, the samples were mounted in clearing solution (50% glycerol, 20% lactic acid diluted in water) and observed by differential interference contrast microscopy using a Zeiss Axioimager Z2 microscope. Picture acquisition was done with a color sCMOS camera (Axiocam 506 color Zeiss).
Statistical analysis
Request a detailed protocolTo identify the main cellular descriptors – cell area, cell volume, sphericity, ellipticity oblate, ellipticity prolate – explaining variance of each cell types, at each ovule primordium growth stages and between genotypes, we used principal component analysis (PCA) on a cells' descriptors matrix. PCA was executed using the R software version 3.6.3. In all different subsets of data, PCA was performed by singular value composition of the centred and scaled data matrix. Data entries with missing values were removed before analysis. All Cell data as exported by Imaris were uploaded but only cells with a ‘Cell Label’ as described in the method section ‘3D imaging’ were plotted. For visualization of PCAs, only the first two principal components were represented in both score and loading plots.
To determine if the means between two datasets were significantly different we used twotailed ttest when the data were normally distributed (n > 30). For all datasets with n < 30, we assume that normality was not possible to assess properly (as small samples most often pass normality tests), thus we used nonparametric tests. Wilcoxon SignedRank twotailed test was used for paired quantifications, and Mann Whitney U twotailed test was used for unpaired quantifications. Tests were performed in Excel or in R (wilcox.test function). To compare ovule proportions, we used twotailed Fischer’s Exact test, using https://www.langsrud.com/fisher.html, available online. Variability was assessed using the Standard Error of the mean (SE), indicated in the graphs and/or the supplemental data when applicable. For some datasets, boxplots were used to improve visualization of data distribution. The number of samples and biological replicates for each experiment are indicated in the figure and/or figure legends and/or supplemental data.
Appendix 1
1.1 Mass springbased ovule growth simulations at cellular resolution
1.1.1 Template preparation and initialization  Regular grid of cells
Following observation of microscopy images of longitudinal sections of carpels showing the placenta, a regular grid of cells has been created with MorphoGraphX (Barbier de Reuille et al., 2015, https://morphographx.org/, Cell Maker plugin). The grid is composed of 28 columns by 5 rows of square cells of side 5.44 μm, the cell rows are staggered with respect to each other by 1.36 μm in the longitudinal direction.
The mass spring mesh is constituted by interconnected segments outlining the cells boundaries and representing the cell walls. The segments have a varying length of 1.36 μm and 2.72 μm, the top margin of the grid of cells has a segment length of 5.44 μm (this is to increase mass spring stability w.r.t. compressive forces generated during ovule growth). A dual mesh connecting the cell centers is associated to the massspring mesh.
The template initialization as well as the proper growth simulations are run with an inhouse built plugin for MorphoMechanX (https://morphographx.org/morphomechanx/) available at a git repository: github.com/GabriellaMosca/MassSpring_2DovuleGrowthModel (Mosca, 2021b).
Once the template has been generated, it needs to be initialized with the parameters illustrated in Figure 3—figure supplement 1B. These properties are assigned at the cell level (no gradient for these properties within a cell), so on the dual graph.
The growth signal is the result of a discrete, cellbased diffusive process. The user sets some cells with fixed high concentration and some other cells with fixed lower concentration of the diffusive substance, while the other cells get their concentration as a result of the diffusive process. More specifically, for the simulations of the growth signal in this work, some cells were assigned a fixed concentration of 0.2 (the vertical column made by the three cells starting from L1 with the highest growth signal concentration in the template, see Figure 3—figure supplement 1B) and some cells were assigned a fixed concentration of 0 (a wellshaped arrangement of cells around the ones with a concentration different from zero, as can be seen in Figure 3—figure supplement 1B). In the model, this is done by running the process Model/Cell Ovule Growth/10 Growth Signal/a Set Cell Type/Assign, after setting the variables to the desired values. The units have been adimensionalized.
Prior to pressurization, Dirichlet boundary conditions in the xdirection have been assigned to the lateral cell walls of the template (this means that those walls are not allowed to displace in the xdirection during pressurization). This condition has been assigned to simulate the presence of the surrounding lateral tissue. Dirichlet boundary conditions can be assigned by running the process Model/Cell Ovule Growth/20 Mass Spring Mechanics/a Set Dirichlet after having replaced the 3D null vector in the field ‘Dirichlet’ with a value different from zero for each coordinate to be fixed. This will affect all the selected nodes of the mass spring mesh.
The process Model/Cell Ovule Growth/20 Mass Spring Mechanics/b Set Cell Type Mechanics and Growth/Assign has been used to assign (i) the different growth directionality (see Figure 3—figure supplement 1B); (ii) different strain threshold for strainbased growth for the L1 cells and the inner ones with a growth signal different from 0 (the remaining cells deform elastically, but are not allowed to perform any growth) (see Figure 3—figure supplement 1B); (iii) different target area for cell division (53 μm^{2} for L1 cells, 100 μm^{2} for the SMC precursor, 80 μm^{2} for the L2 cells neighbors of the pSMC, while all the other cells got the global target area for division of 63 μm^{2}, see Figure 3—figure supplement 1B). To correctly use the process, it is necessary to label L1 cells as ‘L1’ and the other cells with special assigned properties as ‘Special Cell’.
Furthermore, the L1 is prescribed to perform only periclinal growth, by running the process Model/Cell Ovule Growth/20 Mass Spring Mechanics/d Only Periclinal Growth in L1 with the field ‘Value’ set to ‘True’. When this option is assigned, the L1 cell edges corresponding to anticlinal walls are prevented from growing (this process will act on the cells previously labeled as ‘L1’).
It is now necessary to set the polarization field used to define growth directionality: this is obtained as the gradient of a diffusion process. The user needs to select cells with fixed high and low concentration (the specific value assigned here is arbitrary and functional only to be able to set a gradient field which will be normalized). In our simulation the cells in the left side of the template have been set with fixed high concentration, while the cells in the right side of the template have been set with fixed low concentration (whether high/low concentration is fixed on the left or right corner is irrelevant): this produced the gradient visible in Figure 3—figure supplement 1B. To set the cells with fix concentration for the polarization process, the user needs to run the process Model/Cell Ovule Growth/30 Growth Process/a Set Fixed Concentr Diffusion Polarizer/Assign, after having assigned in the dialog box the desired values.
The remaining field values are global and defined in the Table of Figure 3—figure supplement 1B.
After this preparation, the template can be inflated to the prescribed pressure (in our model 0.5 Kg/s^{2}). It is important to stress out that in these simulations, the pressure adopted does not have the units of a physical pressure, but those of a force acting on a 1D segment (instead of on a 2D surface). Connected to this point, we stress out that these 2D simulations have a purely qualitative purpose and no quantitative comparison between the parameters used and the biophysical quantities should be attempted.
After pressurization of the template, new Dirichlet boundary conditions are assigned: the nodes outlining the boundary of the template (except for the L1 top nodes) are blocked in all degrees of freedom (XYZ). This mesh is saved as starting point for the growth and cell division simulation with mass springs.
1.1.2 Confocal microscopybased template preparation
A confocal stack from an almost flat placenta (dataset 2404_275a) has been used. The raw stack has been loaded in MorphoGraphX and by means of the clipping planes, a longitudinal central section of it has been selected. Then, the same procedure described in the MorphoGraphX manual (available here: https://www.mpipz.mpg.de/4085950/MGXUserManual.pdf) for generating $2\frac{1}{2}\mathrm{D}$ mesh has been followed, with the only difference that the surface mesh is generated with the CellMaker plugin as a rectangle with side length 0.5 μm. The signal projection depth has been calibrated so that only the walls of the midsection get projected (the longitudinal mid section has been as well selected accordingly to not display spurious walls from cells at different depth) and set to the value of 1 μ. The final mesh gets converted into a 2D Cell Mesh by the process: Process/Mesh/Cell Mesh/Convert to 2D Cell Tissue with the field ‘Max Wall Length’ set to 10. After cleaning spurious boundary cells, and some manual curation to remove edges which are not junctions on the top wall of L1 cells (this to reduce the artifact of edge collapsing due to the absence of bending resistance during compressive forces), one obtains the final template to be used in the simulation within MorphoMechanX. The template preparation is then identical to what explained in the previous section and the distribution of cells with fixed high concentration is shown in Figure 3—figure supplement 1F.
1.1.3 Simulation cycle
The growth and division simulation starts from the pressurized template at mechanical equilibrium (see following). A simulation cycle consists of the following processes:
growth signal diffusion calculation;
polarization field calculation;
mass spring mechanical equilibrium calculation;
growth step;
cell division.
The simulation cycle is repeated as many times as indicated by the simulation time T (which is not in physical units) mentioned in the main text, in Figure 3 and Figure 3—figure supplement 1. We explain now in detail how each process is computed.
Growth signal diffusion calculation
The growth signal is the result of a diffusive process acting on the cell level (on the dual graph). In this formulation, no sources or sinks are present: only some cells with a fixed concentration are specified at the beginning of the simulation. The discretized diffusion equation, following what done in Bayer et al., 2009, models the concentration variation inside a cell according to this formula (expressed in adimensionalized units through rescaling of the diffusion constant D):
where ${\rho}_{i}$ indicates the signal concentration in the cell $i$, τ, $D$ the diffusion coefficient, ${A}_{i}$ the area of cell $i$, ${N}_{i}$ is the set or cells connected to the cell $i$, while ${l}_{ij}$ is the length of the interface between the two cells. The rationale for this heuristic formula, formally identical to the discrete LaplaceBeltrami operator if the weight factors are made to coincide with ${l}_{ij}$ (Reuter et al., 2009), is that the variation of concentration in the cell $i$ depends on the difference of concentration between this cell and the neighboring ones rescaled by the cell surface (the bigger the cell, the less the concentration variation) and the length of the interface with each neighboring cell. Since we are looking for the steady state solution in the absence of sources and sinks, the rescaled diffusion coefficient $D$ does not affect the final solution and so it is set to one in the following calculations. To find the steadystate concentration, an adaptive forward Euler scheme was adopted (Teukolsky et al., 1992); the problem is easily solvable by a direct method, but this diffusion process has been designed to treat more general, non linear problems as well. So an iteration of the forward Euler scheme, to calculate the variation in concentration for the cell $i$ at the discrete simulation time step $n+1$, is as follows:
The known fixed concentration are initialized with their assigned value. The time increment $\delta \tau $ is chosen based on an adaptive scheme which evaluates the magnitude of the concentration variation at the previous step to determine the new time increment (constrained between a preset max and min allowed time increment). The simulation is considered converged when the update on the concentration ${\rho}_{i}({\tau}_{n+1}){\rho}_{i}({\tau}_{n})$ is smaller than a prescribed tolerance, more precisely when the mean between the averaging of the concentration update over all vertexes and the maximal of the concentration is below a threshold (concentrations taken as absolute values).
Polarization field calculation
The polarization field is computed as the gradient of a diffusive process which is regulated by cells with a fixed concentration prescribed. The evolution process itself is computed as in the previous paragraph. Once the equilibrium concentration has been obtained, its 2D gradient (on the cell grid) is computed and then normalized. This gradient indicates the direction of the ”polarizer’ ${\mathbf{\mathbf{n}}}_{Pol}$, which is used to control the direction of anisotropic growth. The advantage with respect to using fixed coordinates is that the polarization field is updated while the mesh changes, creating a system of local coordinates that follows the organ natural curvature while it grows (Kennaway et al., 2011).
Mass spring mechanical equilibrium calculation
The massspring system is constituted by the network of linear springs which describe the cell walls in the template. The mechanical equilibrium is obtained when the internal reaction forces generated by the mass springs counterbalance the externally applied loads at each point of the mass spring network (MSN):
where ${\mathbf{\mathbf{F}}}_{{I}_{v}}$ indicates the internal reaction force at the node $v$, ${\mathbf{\mathbf{F}}}_{{E}_{v}}$ the sum of the external forces applied to the node $v$ and $v$ belongs to the massspringnetwork ($MSN$). More specifically, the internal forces are given by the Hooke’s law for masssprings:
$k$ is the spring constant (Kg/s^{2}), ${N}_{v}$ is the set of points which are connected to $v$, ${l}_{u\to v}$ is the rest length of spring connecting the point $u$ to the point $v$, while ${\mathbf{\mathbf{p}}}_{u}$ and ${\mathbf{\mathbf{p}}}_{v}$ are the respective positions. Note that $\frac{{\mathbf{\mathbf{p}}}_{u}{\mathbf{\mathbf{p}}}_{v}}{\parallel {\mathbf{\mathbf{p}}}_{u}{\mathbf{\mathbf{p}}}_{v}\parallel}$ is the unit normal vector ${\mathbf{\mathbf{n}}}_{u\to v}$ in the direction connecting the point $u$ to the point $v$. In our case forces are measured in μN and lengths in μm, but, being this a 2D representation of a 3D problem by means of 1D structures such as the mass springs, no quantitative comparison with biomechanical properties of the cell tissue should be attempted.
Regarding the externally applied forces, turgor pressure is considered. In our 2D massspring model this acts on the 1D line segment connecting two vertices in the mass spring system:
here ${\mathbf{\mathbf{F}}}_{{P}_{v}}$ stands for the force due to turgor pressure on the node $v$, ${\mathbf{\mathbf{n}}}_{u\to v}^{\perp}$ indicates the unit vector orthogonal to ${\mathbf{\mathbf{n}}}_{u\to v}$ and pointing outwards with respect to the cell centroid. The summation is indented over all the vertexes $u$ connected to $v$ and over all the cells $C$ whose the edge $u\to v$ is part of the boundary (${\mathcal{B}}_{C}$) and ${P}_{C}$ is the pressure afferent to cell $C$. To compute the forces equilibrium a semiimplicit Euler method has been adopted (Teukolsky et al., 1992):
Equation 1.6 is in a matrix form: $\mathbf{\mathbf{p}}(n)=({\mathbf{\mathbf{p}}}_{1}(n),\mathrm{\cdots},{\mathbf{\mathbf{p}}}_{m}(n))$, ${\mathbf{\mathbf{F}}}_{TOT}(n)=({\mathbf{\mathbf{F}}}_{TOT1}(n),\mathrm{\cdots},{\mathbf{\mathbf{F}}}_{TOTm}(n))$ and $\mathbf{\mathbf{J}}(n)$ is the Jacobian (here computed numerically) of the force vector with respect to the nodal position vector $\mathbf{p}(n)$. $\mathbf{p}}_{i$ is the position of the $i$th node in the $MSN$, while $n$ is the $n$th iteration of the solver. The time stepping, $h$, is computed adaptively and based on the norm of the incremental displacement. The mechanical equilibrium is considered reached when the mean between the averaging of the norm of forces over all vertexes and the maximal force norm is below a threshold. The nodes with prescribed Dirichlet boundary condition will have the corresponding force equal to zero in the direction in which the condition is applied.
In case the option ‘Wall Stiffness based on morphogen’ in the process Model/Cell Ovule Growth/20 Mass Spring Mechanics/c Set Wall Stiffness/Assign has been set to true, for the selected cells during the simulation runtime, stiffness will be assigned within the range specified by the user in a linear negative correlation with respect to the growth signal.
Growth step
After mechanical equilibrium has been reached, an incremental growth step, acting on the rest lengths of the mass springs (${l}_{u\to v}$) is computed. Different types of growth can be specified and added together. In the model there are signal based growth and strain based growth, both can be explicitly required to occur along or orthogonal w.r.t. the polarization field direction ${\mathbf{\mathbf{n}}}_{Pol}$. The incremental growth operator acting on the rest length of a mass spring segment connecting $u$ to $v$ is given as follows:
Here ${G}_{u\to v}^{S}$ is the signal based growth contribution, ${G}_{u\to v}^{\u03f5}$ is the strain based growth contribution, ${G}_{u\to v}^{{\u03f5}_{Pol}}$ is the strain based growth contribution in the polarizer basis ($\mathbf{n}}_{Pol},{\mathbf{n}}_{Pol}^{orth$), $T$ is the discrete simulation time (number of simulation loop iterations) and $\delta \tau $ is a coefficient interpretable as a global growth scaling factor. More specifically:
where ${K}_{Par}^{S}$ is the growth coefficient in the direction parallel to the polarization field, while ${K}_{Per}^{S}$ is the growth coefficient in the orthogonal direction. Since growth is performed by updating a rest length, the growth operator perse can not affect directionality and introduce distortion. Therefore the notion of anisotropic growth is implemented by multiplying a rest length ${l}_{u\to v}(T)$ by ${K}_{Par}^{S}$ if its corresponding vector in the current configuration ${\mathbf{\mathbf{n}}}_{u\to v}$ has a normalized projection along ${\mathbf{\mathbf{n}}}_{Pol}$ bigger than 0.6. Viceversa, if the projection is less than 0.4, the rest length will be multiplied by ${K}_{Per}^{S}$. For the remaining cases, an average of the two coefficients will be adopted. Regarding strainbased growth we have:
where $\u03f5(T)=\left(\parallel {\mathbf{\mathbf{p}}}_{u}(T){\mathbf{\mathbf{p}}}_{v}(T)\parallel {l}_{u\to v}(T)\right)/{l}_{u\to v}(T)$ is the linear strain computed at mechanical equilibrium at simulation loop time $T$ and ${\u03f5}_{thres}$ is the strain threshold. The rest length is not allowed to get bigger than the spring length in the current configuration.
The polarized strain based growth is prescribed as follows:
The growth coefficients used in these formulas are specified cellbased, while growth is performed edgebased: therefore the coefficients actually used are an average quantity over the cells related to the edge.
Cell division
Cell division is performed according to the rule ‘shortest wall through the centroid’ as described in Mosca et al., 2018. The target cell areas for division are assigned as described in Sec. 1.1.1. Cell wall pinching has been globally set to 0.1 (which corresponds to a proportionality factor connected to the new wall segment length or the distance between the newly inserted point and the cell centroid, depending on which is smaller), and is constrained to be overall smaller than 0.1 μm. The minimal distance from a preexisting wall has been set to 0.8 μm. Cell wall sampling to search for the shortest wall has been set to 0.05 μm (the starting point for the sampling is aleatoric and in principle not guaranteed to be the same at different runs).
Upon cell division, cell properties (such as fixed growth signal concentration, fixed polarizer, cell type, target division area, special pressure, growth coefficients, etc.) are inherited by the daughter cells. For what concerns fixed high growth signal concentration, selective rules have been implemented to allow it to be propagated along a vertical line in the inner tissue. For what concerns cells in the inner layers, if a dividing wall is in the anticlinal direction, only one of the two daughter cells will inherit the fixed high signal concentration (the one with the smallest distance from the closest cell with a fixed high concentration). If the newly inserted wall is in the periclinal direction, if the mother cell (prior to division) is between two cells with a fixed high concentration, both daughter cells will inherit this property, otherwise only the one connected to the vertical stream of cells with a fixed high concentration. Regarding L1 cells, the daughter cells inherit both the fixed high concentration if the mother cell has two L1 neighbors with a fixed high concentration, otherwise only one cell will inherit this property with preference for the daughter cell with the shortest distance from the L2 cell with fixed high signal concentration. The propagation rules for fixed high concentration are not equivalent to a proper canalization model (which was not the aim of these simulations) and sometimes manual curation during the simulation time evolution has been required to ensure an uninterrupted vertical line of cells with fixed high concentration.
To assess the stability of the reference simulation (MSModel 2), a robustness analysis has been performed by varying independently stiffness and growth signal intensity of 10% around the reference value. The results are qualitatively compatible with the reference model and available in the git repository: github.com/GabriellaMosca/MassSpring_2DovuleGrowthModel/tree/master/robustnessAnalysis.
1.2 Cell anisotropy quantification
Similarly to what done in Louveaux et al., 2016, a shape matrix has been defined for each cell. This matrix is the covariance matrix of the coordinates of the points along the cell outline. Given a distribution of points in the 2D plane, their covariance matrix with respect to their coordinates in the Cartesian plane is defined as
where $x}_{i$, $y}_{i$ are the ith point coordinates, ${N}_{c}$ is the number of points and x_{0}, y_{0} are the arithmetic average of all the points coordinates.
This matrix is symmetric and so diagonalizable, its eigenvectors indicate the directions of maximal and minimal points dispersion (and so the anisotropicity of the distribution). More precisely, the eigenvector connected to the biggest eigenvalue defines the direction around which the standard deviation of the projection of the points position along that axis is maximal, while the other eigenvector defines the direction around which this same quantity is minimal. At the same time, one can see the first eigenvector as the one minimizing the variance of the distance of the points distribution from the direction it indicates, while the second eigenvector as the one maximizing instead this quantity.
If the points are uniformly distributed around a shape outline, the eigenvectors can be used as a proxy for the maximal and minimal elongation directions of the shape, in the sense described above. If one considers a rectangle (where the concept of shape anisotropy can intuitively be assimilated to the aspect ratio of its sides), the eigenvectors of this covariance matrix will be parallel to its sides and, if the sides ratio is greater than 1, also the corresponding eigenvalues will be greater than one. Furthermore the eigenvalues ratio depends on the sides lengths, but not on the rectangle area. If a rectangle has an aspect ratio greater than another (always computing longest size over shortest side), this relation will be valid also for the corresponding eigenvalues ratio.
If a rectangle is now a 2D cell which is growing and it is getting more elongated while not changing in the other direction, the eigenvalues ratio will increase. On the other hand if, while being anisotropic, it is still growing in an isotropic manner, its eigenvalues ratio will not change. So by looking at the eigenvalues ratio and their variation over time, it is possible to deduce the cell anisotropicity regardless of its surface and how it varies over time.
When a shape differs from that of a rectangle (or an ellipse), it is not possible to use the sides (main axes) aspect ratio as a measure of anisotropy, but the shape matrix provides an extension of this concept to other shapes. Clearly, the more the shape deviates form that of a rectangle, the more care should be used to interpret the results (especially with curved shapes, where a nonlinear approach should be adopted).
In the measurements performed in this work, the cell perimeter has been enriched with intermediate points at an average distance prescribed by the user and the averaging has been weighted with the points distance to avoid artifacts due to non perfect homogenous points distribution. So the actual equation used to compute the covariance matrix is:
Where ${P}_{c}$ stands for the cell perimeter and $l}_{i$ for the length of the perimeter segment containing the node $i$ in its middle. Furthermore, for the images shown in this work, it has been selected the option to measure cell anisotropy along the two orthogonal growth directions of the ovule. For this reason, by using ${\mathbf{\mathbf{n}}}_{Pol}$ (the vector parallel to the polarization field direction, see Sec. 1.1.3) and ${\mathbf{\mathbf{n}}}_{Pol}^{\perp}$, the orthogonal one, these two quantities are computed:
and their ratio is used as a proxy of cell anisotropy along the organ growth direction and the one orthogonal to it. More precisely, as in the figures the aim was to show the cell elongation along the ovule main growth direction, the anisotropy ratio displayed is given by this formula: ${\alpha}_{Pol}^{\perp}/{\alpha}_{Pol}$. This means that a displayed value of one indicates that the cell is fully isotropic, a displayed value bigger than one indicates that the cell is more elongated along the organ main elongation axis (${\alpha}_{Pol}^{\perp}$), while a value between one and zero indicates that the cell is more elongated in the orthogonal direction to the main organ elongation axis. This calculation is performed in the model by running the process: Model/Cell Ovule Growth/50 Shape Quantifier/01 Compute Skew Symmetric Tensor, where the user can specify the average distance for the virtual points to be placed along the shape perimeter for the computation (one should refine this number until the result is not significantly affected). The eigenvectors displayed in Figure 3E are scaled by their pertinent eigenvalue divided over the sum of all three eigenvalues (normalization), all globally rescaled by a length factor.
1.3 2D FEM based simulations for ovule growth  continuous limit
The finite element method (FEM) simulations are handled in MorphoMechanX by an inhouse developed tool available at the git repository github.com/GabriellaMosca/FEM_2DOvuleGrowthModel (Mosca, 2021c). They run on a 2D mesh tessellated with triangles which constitute linear membrane elements with a mathematically prescribed thickness (not relevant in this type of simulation though).
1.3.1 Template preparation
The initial template is build as rectangle tessellated by squares assembled in 28 rows and 56 columns. The squares have a side of length 13.0527 μm (the value assigned to the length is arbitrary as these simulations aim at a qualitative analysis in 2D of a 3D problem). The tessellation squares have been each triangulated with four triangles sharing a vertex in the square center (see Figure 3—figure supplement 1A).
This template represents a portion of the ovule placenta as seen through a longitudinal cut. The first three lines or squares from the top are identified with the L1 layer and are assigned specific properties different from the remaining tissue, as will be described in the following.
The description which follows is referred to the reference case template, FEMmodel 1, as from Table 1 in the main text. The other models are obtained through minor modifications of the process described in the following.
As first thing, the growth signal intensity is assigned: this is obtained, similarly to what done with the mass springs, as the outcome of a diffusive process with some elements at a high and low fixed concentration. First of all, all the faces of the mesh need to be selected and the process Model/CCF/Fem Diffusion Growth/04 Create Diffusion Element has to be run, so that the element for growth signal diffusion is generated.
Then the nodes on the mesh with fixed signal concentration (high and low) need to be set. The nodes with high fixed growth signal concentration are forming a vertical line made of 16 rows by seven columns starting from the top margin of the mesh, in a centered position (see Figure 3—figure supplement 1A for the signal distribution). By running the process Model/CCF/Fem Diffusion Growth/03 Set Diffusion Dirichlet with the field ‘Value’ set to 0.2 (adimensional units adopted for diffusion as in the case of mass springs), the fixed high signal concentration is assigned to these vertices. Note that the ‘Dirichlet Attribute’ and ‘Morphogen Attribute’ (as well as the ‘Element Attribute’) names need to be assigned consistently for all the subprocesses contained in Model/CCF/Fem Diffusion Growth.
The vertices with a fixed low signal concentration are assigned with the same procedure. These vertices form the perimeter of a rectangle which encloses the vertical strip of cells with fixed high concentration. It is located at a distance from the mesh middle (in the horizontal coordinate) of 8 vertices, and from the mesh top (in the vertical coordinate) of 21 vertices (see Figure 3—figure supplement 1A for the distribution of the signal). The fixed assigned concentration is 0.
After assigning the initial conditions, the diffusion process itself is run by Model/CCF/Fem Diffusion Growth/00 Diffusion Solver. The details on how the diffusive process is handled are provided in the paragraph 1.3.2.
Similarly to what done for the growth signal, the material and growth main directions (E2, KPar) are obtained through a diffusive process (with the exception of the L1 layer, which is set with a different method). More precisely the gradient of the diffusive process will be used, or, alternatively, the direction orthogonal to it. The process for assigning the directions is Model/CCF/Fem Diffusion Anisotropy and is formally identical to Model/CCF/Fem Diffusion Growth.
It is possible to handle growth and material anisotropy with a unique diffusive process, where the source with fixed high signal concentration is located at the central vertex (in the horizontal direction) at the mesh top, while the vertexes with fixed low concentration are located at the mesh boundary (excluding the top boundary): the outcome of this diffusive process produces a radial gradient pointing in the direction of the vertex with fixed high concentration (see Figure 3—figure supplement 1A, where the field is visible for the inner layers). For the mesh portion corresponding to the L1 layer, the material and growth anisotropy are set with a separate process, as it will be explained in the following.
At this point comes strictly mechanical part of the simulation. The mechanical elements need to be created by selecting all the faces of the mesh and by running the process Model/CCF/04 Reference Configuration which not only generates the elements for the FEM solver, but also sets their reference configuration as the current one. This process assigns also a thickness value, but this parameter is irrelevant in this 2D membrane simulation (so its has been left to 1 μm, as the membrane thickness won’t affect the equilibrium configuration or the stresses, but will create a global rescaling factor for resultant forces).
Subsequently material properties get assigned. The template has a linear, hyperelastic, transverse isotropic, St. Venant material law, characterised by a special fiber direction (in this text named E2); for a detailed description of this material law and how it is used to compute the force vector for mechanical equilibrium (Hofhuis et al., 2016).
As first thing, it is necessary to assign the intensity of the Young modulus in the fiber direction (E2) and in the isotropic plane (denominated E1E3 in the text: in the template this is represented by a vector, but conceptually it is a plane because of the virtual thickness). The process Model/CCF/08a Set Aniso Dir sets global values for the Young modulus and Poisson ratio. For the described template, it needs to be run with the field ‘Young E1E3’ equal to ‘100’, ‘E2’ to ‘500’ and ‘Poisson’ equal to ‘0.3’. An additional stiffness, negatively proportional to the growth signal previously computed through diffusion, is assigned to the E1E3 plane by running Model/CCF/06b Material Properties based on Morphogens (for faces) in all the models as it creates a sharper delimitation of the ovule protrusion w.r.t to the placenta (for a comparison see the outcomes at the git repository github.com/GabriellaMosca/FEM_2DOvuleGrowthModel/tree/master/constantStiffness)
where ${Y}_{E1E3}$ is the Young modulus in the isotropic plane, γ stands for the proportionality factor (set to 400), and growth signal is the signal obtained from the diffusive process Model/CCF/Fem Diffusion Growth(make sure to put the proper field signal name in the dialog box, its name was specified in Model/CCF/Fem Diffusion Anisotropy/05 Morphogen Visualize). The Young modulus distribution is the one shown in Figure 3—figure supplement 1A. In the mass springs simulation the stiffness was instead assigned constant in the template, as a replication of what done with the FEM simulation produced a numerically unstable model (see github.com/GabriellaMosca/MassSpring_2DovuleGrowthModel/tree/master/variableStiffness). We argue that for mass spring models, the effect of softer cell wall is already implicitly encoded in the combination of signal based growth and passive strain based growth, so that adding such an effect explicitly exaggerates the dishomogeneity between faster and more slowly growing cells. In the FEM simulation, lacking cellular resolution, intead, the passive strain based growth is mostly acting as a relaxation to accumulation of stresses due to tissue conflicts.
For the all the template, except for the L1 layer, the special fiber direction (E2) is given by the vector orthogonal to the gradient of the diffusive field computed in Model/CCF/Fem Diffusion Anisotropy. To assign it this way, after selecting the face elements of the template, the user needs to run Model/CCF/08b Set Aniso Dir Morphogens, where, in the field ‘Direction Type’ the key ‘E2’ has been assigned, in the field ‘Direction from Diffusion Gradient’ the key ‘Orthogonal to Gradient’ and in the field ‘Morphogen Signal’ the name chosen for the morphogenetic signal in Model/CCF/Fem Diffusion Anisotropy.
The L1 layer gets assigned a global E2 direction coaligned with the yaxis (the vertical coordinate) by running the process Model/CCF/08a Set Aniso Dir with the accuracy of specifying direction type as ‘E2’ and ‘Direction’ as ‘0 1 0’.
The same procedure is repeated to assign the growth anisotropy direction (Mosca et al., 2018), denominated KPar in this text, with the difference that the field ‘Direction Type’ is ‘KPAR’ and the ‘Direction from Diffusion Gradient’ is ‘Parallel to Gradient’. Regarding L1 layer, this is assigned a ‘Direction’ of ‘1 0 0’.
An edgepressure acting on the template boundaries has been assigned (as the lateral and bottom boundaries will be constrained by Dirichlet condition, it is actually assigned only to its top boundary). This edge pressure acts exactly like the one described for the mass springs in Equation 1.5 and has the same units. It has been set to a value of 5 Kg/s^{2} and its role is to simulate, to some extent, the action of turgor pressure into generating stresses in the tissue.
Dirichlet boundary conditions (as with mass springs) are assigned on the bottom and lateral boundary of the template by running the process Model/CCF/15 Set Dirichlet and blocking all the degrees of freedom (setting the field ‘Dirichlet Label’ to ‘1 1 1’).
As last, it is required to assign the growth field. In the reference simulation (FEMmodel 1) signal based growth is purely in the polariser direction (KPar) combined with strainbased growth (see following for explanation). This is assigned for the selected faces by running the process Model/CCF/11b Set Growth Morphogens (on Faces) after setting the field ‘KPar Scaling Factor’ to 1, all the others to 1, and properly assigning the name in the field 'Morphogen Signal 1’.
1.3.2 Diffusion simulation
Diffusion simulation is used in this model to generate a concentration field given some fixed concentration nodes. The problem is solved in 2D on a flat surface. It is slightly different from the one implemented for mass springs (see paragraph 1.1.3) as it is based on the FEM and starts from a continuous formulation originally expressed by this equation (heat equation):
where ∇ is the Laplace operator, $c$ the concentration of the diffusible substance and $\mathbf{\mathbf{x}}$ are the coordinates used to described the space where diffusion is occurring. This problem requires boundary conditions (otherwise its solution would be trivial). For this diffusive process, unless prescribed Dirichlet condition (fixed concentration of the diffusing substance) have been assigned, zero flux condition is naturally assumed along the boundary. Internally to the domain, fixed concentration can be assigned as well, even if mathematically, this is equivalent to create an internal boundary (i.e. if inside a connected 2D surface a circle gets prescribed fixed concentration, which does not need to be constant in space, what will affect the solution outside this circle is only its boundary and not what has been set inside).
An explanation or introduction to the FEM is outside the purpose of these supplementary information, but there are several sources and manuals to consult, as for example Hughes, 2012. In this case a direct solver method has been adopted, which means that the problem can be written, after the FEM discretisation, in the following way:
where $\mathbf{\mathbf{J}}$ is a matrix, $\mathbf{\mathbf{c}}$ is the vector of the nodal concentrations and $\mathbf{\mathbf{B}}$ is a vector. The matrix $\mathbf{\mathbf{J}}$, is given, in components, is as follows:
where the external sum is over all the triangles $\U0001d517$ of the FEM mesh containing the node $i$, while the internal sum is over all the nodes $j$ which are connected to the node $i$. $A}_{\mathfrak{T}$ is the triangle area, $\mathbf{\mathbf{D}}{\varphi}_{i}$ is the vector of derivatives (in 2D) of the triangle linear basis functions ${\varphi}_{i}$ (Hughes, 2012) referred to the triangle $\U0001d517$ (which has been omitted for simplicity of notation) and node $i$, explicitly:
with $X$ and $Y$ the local planar coordinates of the triangle. Regarding the B vector, it can be computed in components as follows:
the final nodal concentration vector $\mathbf{\mathbf{c}}=({c}_{1},\mathrm{\cdots},{c}_{m})$ is obtained as the solution to the matrix equation $\mathbf{\mathbf{J}\mathbf{c}}=\mathbf{\mathbf{B}}$ after adding to it the fixed nodal concentration values.
In this calculation the diffusion coefficient $D$ has been omitted as it does not affect the steady state concentration in this type of diffusion.
1.3.3 One simulation cycle
A simulation cycle, which gets repeated in loop as many times (T) as specified by the user, consists of the following units:
mechanical equilibrium simulation;
growth step;
mesh subdivision.
The template gets altered from a mechanical equilibrium condition because of the action of the edge pressure (assigned to the top margin) and of growth.
Mechanical equilibrium simulation
After pressurization and each growth step, the system generally needs to find a new equilibrium configuration which accommodates the pressure induced stresses as well as the residual ones introduced by the growth process (Rodriguez et al., 1994; Goriely and Ben Amar, 2007 for an explanation of residual stresses).
As already introduced in Sec. 1.3.1, the material is transversely isotropic with a SainVenant Kirchhoff material law.
The equilibrium is computed by minimizing its strain energy by means of a semiimplicit Euler scheme (see Equation 1.6). The elemental contribution to the nodalforce vector calculation is provided in detail in Hofhuis et al., 2016. The final nodalforce vector is obtained by summing for each node all the elemental contributions it is connected to. The Jacobian of the force vector is obtained through numerical derivation.
Details of material parameters and numerical convergence are provided in the Table in Figure 3—figure supplement 1A.
Growth step
Growth is performed as a series of small deformations which are applied locally to the reference configuration and deform it without, directly, inducing stress or strain. As the deformation is local, no compatibility requirement is naturally satisfied and this has to be enforced by combining the growth deformation (which does not cause strain or stress) with an elastic deformation (which might induce strain and stress even in the absence of external loads), see Rodriguez et al., 1994.
The underlying hypothesis in this approach is that the timescale of growth is much bigger than that of elastic waves propagation, so that growth is supposed to occur in a quasistatic condition of mechanical equilibrium (Goriely and Ben Amar, 2007).
In this formulation it is assumed (with no loss of generality) that all the rotational components of the growth tensor are absorbed in the elasticity tensor, so that growth is described purely by a stretch tensor, with no rotational component. In this description of growth, for each growth step, the reference configuration is the one obtained at the previous growth step and not the original one and it is called virtual as it might not ever occur physically (because of incompatibility). So, if ${\mathbf{\mathbf{F}}}_{E}^{n}$ and ${\mathbf{\mathbf{F}}}_{G}^{n}$ are the elastic deformation tensor and the growth deformation tensor at the $n$th iteration, if ${\mathbf{\mathbf{X}}}_{n1}$ is the virtual reference configuration obtained at growth step $n1$, the global deformation ${\mathbf{\mathbf{F}}}_{EG}$ at the $n$th growth step will be as follows:
In general the growth tensor ${\mathbf{\mathbf{F}}}_{G}^{n}$ might depend on the current body coordinates ${\mathbf{\mathbf{x}}}_{n1}$ or on its state of elastic deformation/stress as computed from the previous virtual reference configuration. In the simulations in this paper, growth is constituted by an anisotropic part explicitly specified by the user through two diffusive process (one for the anisotropy direction ${\mathbf{\mathbf{v}}}_{Kpar}$ and one for the growth intensity $KPar$, as described in Sec. 1.3.1) and a strainbased growth.
Growth is performed locally, by updating the rest length of the elemental triangles of the reference (possibly incompatible) configuration. It is possible to update this approach because triangles are simplexes, so their shape is fully determined by their sidelengths. The physical position of the reference triangle in space is not relevant.
Each triangle in the reference configuration stores in a ordered way the rest lengths of its edges and the angles made by the material fibers and growth anisotropy direction with its first side (the sides ordering is consistently preserved during the simulation). At the very first growth step, the growth direction and material anisotropy direction are naturally expressed in the reference configuration, so it is just a matter of computing the angle they make locally with the first triangle side. The growth anisotropy angle will be denominated ${\alpha}_{KPar}$, while the material anisotropy angle ${\alpha}_{m}$. The reference triangle is initially built from its ordered rest lengths in a local frame, where its first side is coaligned and positively oriented with the xaxis. In this way
The growth tensor for this growth component, ${\mathbf{\mathbf{G}}}_{KPar}$, is expressed in a diagonal basis:
where ${g}_{Kpar}$ is the growth coefficients for the local anisotropy growth in the direction parallel to ${\mathbf{\mathbf{v}}}_{Kpar}$. In general it is possible to specify growth also in the direction orthogonal to ${v}_{KPar}$, but, depending on the growth hypotheses made, this direction might be orthogonal to ${v}_{KPar}$ in the current configuration, but not in the reference one. Or it might lose the notion of orthogonality in both configurations if it is assumed that growth fields are specified at the beginning of the simulations and get deformed during growth and elastic deformation. This of course, only in the case anisotropic growth is coupled with another growth tensor which makes the final growth operator nondiagonal w.r.t to ${v}_{KPar}$ anymore. It will not be exposed here how the growth tensor should be computed in these cases.
Consequently, it is necessary to rotate the reference triangle so that ${\mathbf{\mathbf{v}}}_{Kpar}$ is representable in a canonical basis coordinates: given the explicit form of ${\mathbf{\mathbf{G}}}_{KPar}$, as (0, 1).
The strainbased growth tensor ${\mathbf{\mathbf{G}}}_{E}$, which will be summed up with the purely anisotropic one (${\mathbf{\mathbf{G}}}_{KPar}$) to constitute a global growth tensor ${\mathbf{\mathbf{G}}}_{TOT}$, is expressed in the diagonal basis of ${\mathbf{\mathbf{G}}}_{KPar}$. The strain here considered is the GreenLagrange strain tensor, which is expressed in the current body configuration (more specifically it maps a vector in the reference configuration to a vector in the current one). The computation of ${\mathbf{\mathbf{G}}}_{E}$ components proceeds as follows:
the triangle represented in the current 3D configuration gets mapped by a rotation matrix in 2D;
the GreenLagrange strain tensor between the current triangle (mapped in 2D) and the reference one is computed (see Hofhuis et al., 2016) and eigenvalues and eigenvectors are derived;
after computation of the engineering strain for each inplane eigenvalue, the strainbased growth matrix in diagonal representation is built:
(1.23) ${\mathbf{G}}_{E}^{D}={g}_{\u03f5}\left(\begin{array}{cc}\hfill max(0,{\u03f5}_{11}^{E}{\u03f5}_{0}^{E})\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill max(0,{\u03f5}_{22}^{E}{\u03f5}_{0}^{E})\hfill \end{array}\right)$where ${\u03f5}_{11}^{E}$ and ${\u03f5}_{22}^{E}$ are the ordered engineering strain eigenvalues obtained from the GreenLagrange strain eigenvalues ${\u03f5}_{11}$ and ${\u03f5}_{22}$ (${\u03f5}_{ii}^{E}=1+\sqrt{(1+2{\u03f5}_{ii})}$), while ${\u03f5}_{0}^{E}$ is a threshold value below which no strain based growth occurs and ${g}_{\u03f5}$ is a growth proportionality factor (for each principal direction; the engineering strain has been chosen due to its linear dependence on the stretch);
the strainbased growth tensor gets represented in the same basis which diagonalizes ${\mathbf{\mathbf{G}}}_{KPar}$:
(1.24) ${\mathbf{G}}_{E}={\mathbf{R}}_{eigen}^{T}{\mathbf{G}}_{E}^{D}{\mathbf{R}}_{eigen}$where ${\mathbf{\mathbf{R}}}_{eigen}^{T}$ is the matrix formed by the GreenLagrange strain eigenvectors as columns.
The final growth tensor is then:
and it is expressed in the basis which diagonalises ${\mathbf{\mathbf{v}}}_{Kpar}$. By applying it to the reference triangle coordinates, the triangles is grown (at the end only the rest lengths are stored, as the triangle arrangement in space is irrelevant).
It is now necessary to update the growth anisotropy and material angle with respect to the grown triangle first side, ${\alpha}_{KPar}$ and ${\alpha}_{m}$ respectively. The reference triangle configuration has been rotated so to be in the diagonal basis for ${\mathbf{\mathbf{G}}}_{KPar}$, so the same rotation has to e applied to ${\mathbf{\mathbf{v}}}_{m}$:
where by definition ${\mathbf{\mathbf{R}}}_{KPar}$ is the rotation matrix which rotates ${\mathbf{\mathbf{v}}}_{Kpar}$ into the vector represented by the coordinates (0, 1). At this point, by applying the final growth operator ${\mathbf{\mathbf{G}}}_{TOT}$ to the two vectors, one obtains their orientation (plus an irrelevant stretch) in the new grown configuration, so it is possible to compute the updated angles between these two vectors and the first, consistently oriented, triangle side.
Underlying this approach there is the idea that the anisotropy direction got computed at the beginning of the simulation and gets passively deformed by the growing rest configuration. An alternative approach would have been to recompute the diffusion process and its gradient at each growth step on the current configuration and then map it on the rest configuration. In this specific simulation scenario, it is possible to verify that the two gradients so obtained do not differ significantly.
Mesh subdivision
When a triangle area in the current configuration exceeds a preset threshold (Max Triangle Area in Table in Figure 3—figure supplement 1A), it gets divided by a bisection algorithm which splits its longest side in two segments of equal length and generates so two new triangles (which replace the old one) by connecting the newly inserted vertex to the opposed one. To ensure a conformal mesh, the subdivision is propagated to neighbor triangles as described in Rivara and Inostroza, 1995. After each subdivision the triangle, edge and vertex properties need to be propagated
the vertex properties which need to be propagated to the newly inserted vertex are the Dirichlet condition (both for mechanics and diffusion): if the new vertex has been inserted between two vertexes with Dirichlet conditions assigned, it will inherit them, otherwise no;
an edge in the current configuration gets divided in two equal parts, so the new edges inherit half the rest length of the undivided one each;
by splitting an edge, the initial triangle is split in two smaller triangles (not necessarily identical): each of them inherit the same triangle properties of the undivided one (i.e. material properties, growth properties, anisotropy angles, etc.).
Robustness analysis
As done with mass springs, a robustness analysis around the reference model (FEMModel 1) has been performed by varying independently stiffness and intensity of growth signal of 10% around the reference values. As can be observed by examining the results in the git repository github.com/GabriellaMosca/FEM_2DOvuleGrowthModel/tree/master/robustnessAnalysis such variations do not affect significantly the outcome, confirming the stability of the reference model.
1.4 3D analysis of ovule morphology and connectivity
1.4.1 Ovule segmentation export from IMARIS
A inhouse tool (ExportImarisCells, git repo github.com/barouxlab/ExportImarisCells, Mosca, 2021a) has been developed to export segmentations performed in the proprietary software IMARIS (imaris.oxinst.com). This tool, to be launch within IMARIS, creates a cellmask with a different integer value for each segmented cell and the value zero for everything else. All the cellmasks are summed algebraically together to provide a final mask where each cell (each 3D voxel occupied by a cell) is labeled with a different integer value and the outside space is labeled with the value 0. This algorithm is based on the hypothesis that cells occupy mutually exclusive space regions. It is possible to export with different labels 2^{16} different cells. The output of this operation is saved as a tiff images series.
1.4.2 3D mesh cellbased and surface mesh creation in MorphoGraphX
After exporting the segmentation as a tifffile series from IMARIS in the previous step, it is now possible to load it in MorphoGraphX (Barbier de Reuille et al., 2015, https://morphographx.org/). The user needs to specify manually the z voxel size. At this point the stack is loaded and visible in MorphoGraphX and the option 16 bit as well as label need to be selected to proceed properly. Upon running the process Process/Stack/Segmentation/Relabel the image series will be relabeled and can be saved as a mgxs file, a native format in MorphoGraphX for segmented stacks. After this step, it is possible to properly generate the cell mesh by running the process Process/Mesh/Creation/Marching Cubes 3D, where the cube size, which is related to the mesh refinement, has to be specified by the user (a size of 0.5 ${\mu}^{2}$ has been specified for meshes connected to ovules till Stage 1I, for later stages, to avoid underestimation of surface of contact between cells, a size of 0.1 ${\mu}^{2}$ has been adopted). This provides a cell mesh which can already be used for volume quantifications. As a further step, to be able to use it in MorphoMechanX, the cell mesh needs to be converted and saved as a 3D tissue by running the process: Process/Mesh/Cell Mesh/Convert to 3D Cell Tissue and then saving the result. Both mesh types can be saved in the native MorphographX mesh format mgxm. Starting back from the saved relabeled stack (the mgxs file), by following the instructions provided in the MorphoGraphX user manual (available here https://www.mpipz.mpg.de/4085950/MGXUserManual.pdf) it is possible to generate a surface mesh and by running Process/Mesh/Signal/Project Mesh Curvature the curvature signal will be projected on it. The curvature used in the analysis reported in Figure 1E–F is the minimal curvature, with manually set max and min range (for Stages 0I and 0II the range has been set between −0.15 and 0.15, while for later stages, as the ovule profile gets sharper, between −0.25 and 0.25).
1.4.3 Connectivity analysis for the ovule layers
The 3D Cell Tissue mesh saved as described in the previous section can be loaded in MorphoMechanX to perform a semiautomated cell layer classification and connectivity analysis. An inhouse tool to be run within MorphoMechanX has been developed to be able to detect semiauthomatically cell layers, after the user selects manually the L1 cells. If the user selects as well explicitly the pSMC/MMC, this will be marked as so in the analysis ad the companion cells will be identified (both cell types are considered subcases on L2 layer cell classification). The tool is available at the git repository: github.com/GabriellaMosca/3DAutoLabelingShapeQuant (Mosca, 2021d).
A cell is identified to be a L2 cells if it shares a portion of its wall with the L1 layer. To avoid spurious contacts due to segmentation imperfection, a minimal threshold for contact can be specified (in this case 0.01 ${\mu}^{2}$). Once the L2 cells have been identified (the pSMC/MMC and companion cells are considered L2 cells), the L3 cells will be detected as the cells sharing a portion of their wall with L2, but none with L1. Finally all other cells are gathered in a global group. This tool will save to a file, for each selected cell, its labeling (L1, L2, pSMC/MMC, L3, etc.), its volume and in case of L2 cell (including pSMC/MMC cell) its surface of contact with L3 layer. In the file are furthermore saved the three eigenvalues of the positional covariance matrix described in the next section. For a tutorial on how to use the tool, see the README.md file in the git repository.
1.4.4 MMC shape characterization with its positional covariance matrix
Following the idea presented in Sec. 1.2 of this document, a 3D version of the covariance matrix of Equation 1.11 has been used to characterize anisotropy of pSMC/MMCcell shape in 3D. The computation, analogously with what done for the 2D case, is performed along the cell surface and the already existing mesh is used for the point distribution. The points used for the computation are the centroids of each mesh triangle and the covariance matrix, with elements area as weights to compensate for dishomogeneities in the mesh looks like:
in this formula $N$ is the total element number (not the number of nodes), ${A}_{i}$ is the element area, ${x}_{i},{y}_{i},{z}_{i}$ are the elements centroids and ${x}_{0},{y}_{0},{z}_{0}$ are the overall 3D shape centroid (computed as the weighted average of the element centroids positions).
In this case the three eigenvectors of the covariance matrix correspond to (i)the axis where the variance of the points projection on it is maximal, (ii) the axis where it is minimal and (iii) a third orthogonal direction to the other two where the variance is intermediate.
To get a proxy of cell anisotropy (so independent from the cell size) it is possible to rescale the eigenvalues of the covariance matrix by dividing each for the sum of all three (making them a partition of unity):
where c_{i} are the eigenvalues of the covariance matrix. Furthermore, to be able to compare the principal anisotropy directions of the analyzed cells, the eigenvalues were saved in a rigid Cartesian frame coaligned with the ovule main elongation axis.
Data availability
The data generated or analysed during this study are included in the manuscript and supporting files. Source data files have been provided for all Figures. Segmented Image data are deposited in Dryad Digital Repository, available at: https://doi.org/10.5061/dryad.02v6wwq2c. Cell Statistics from these images are deposited as 'Segmented Dataset.csv' at https://github.com/barouxlab/OvuleViz (copy archived at https://archive.softwareheritage.org/swh:1:rev:fd614aa1e80258928ee036191f26c3dd703d3141).

Dryad Digital RepositoryArabidopsis ovule primordium stages 0I to 2II, wildtype Col_0, Ws4 and botero (bot17) mutant.https://doi.org/10.5061/dryad.02v6wwq2c

GitHubID fd614aa. Segmented_Dataset_Hernandez_Mosca_etal_2020.
References

Megasporogenesis in Arabidopsis thaliana L.: an ultrastructural studySexual Plant Reproduction 12:99–109.https://doi.org/10.1007/s004970050178

Integration of transportbased models for phyllotaxis and midvein formationGenes & Development 23:373–384.https://doi.org/10.1101/gad.497009

Flowers under pressure: ins and outs of turgor regulation in developmentAnnals of Botany 114:1517–1533.https://doi.org/10.1093/aob/mcu187

Bringing Open Data to Whole Slide ImagingDigital Pathology 2019:3–10.https://doi.org/10.1007/9783030239374_1

A computational framework for 3D mechanical modeling of plant morphogenesis with cellular resolutionPLOS Computational Biology 11:e1003950.https://doi.org/10.1371/journal.pcbi.1003950

BookArabidopsis: An Atlas of Morphology and DevelopmentNew York: SpringerVerlag.https://doi.org/10.1007/9781461225980

A continuous growth model for plant tissuePhysical Biology 13:065002.https://doi.org/10.1088/14783975/13/6/065002

Coordination of Morphogenesis and CellFate Specification in DevelopmentCurrent Biology 27:R1024–R1035.https://doi.org/10.1016/j.cub.2017.07.010

Growth and biomechanics of shoot organsJournal of Experimental Botany 70:3573–3585.https://doi.org/10.1093/jxb/erz205

Nuclear envelope: a new frontier in plant mechanosensing?Biophysical Reviews 9:389–403.https://doi.org/10.1007/s1255101703026

Phyllotaxis as geometric canalization during plant developmentDevelopment 147:dev165878.https://doi.org/10.1242/dev.165878

On the definition and modeling of incremental, cumulative, and continuous growth laws in morphoelasticityBiomechanics and Modeling in Mechanobiology 6:289–296.https://doi.org/10.1007/s1023700600657

The molecular and genetic basis of ovule and megagametophyte developmentSeminars in Cell & Developmental Biology 9:227–238.https://doi.org/10.1006/scdb.1997.0214

The developmentalgenetics of canalizationSeminars in Cell & Developmental Biology 88:67–79.https://doi.org/10.1016/j.semcdb.2018.05.019

A mechanical feedback restricts sepal growth and shape in ArabidopsisCurrent Biology 26:1019–1028.https://doi.org/10.1016/j.cub.2016.03.004

Female germ unit: organization, isolation, and functionInternational Review of Cytology 140:233–293.https://doi.org/10.1016/S00747696(08)610992

BookThe Finite Element Method: Linear Static and Dynamic Finite Element AnalysisCourier Corporation.

Molecular regulation of the mitosis/meiosis decision in multicellular organismsCold Spring Harbor Perspectives in Biology 3:a002683.https://doi.org/10.1101/cshperspect.a002683

Molecular mechanisms of robustness in plantsCurrent Opinion in Plant Biology 16:62–69.https://doi.org/10.1016/j.pbi.2012.12.002

Establishing a framework for female germline initiation in the plant ovuleJournal of Experimental Botany 70:2937–2949.https://doi.org/10.1093/jxb/erz212

Katanin: A Sword Cutting Microtubules for Cellular, Developmental, and Physiological PurposesFrontiers in Plant Science 8:1982.https://doi.org/10.3389/fpls.2017.01982

KATANIN 1 Is Essential for Embryogenesis and Seed Formation in ArabidopsisFrontiers in Plant Science 8:728.https://doi.org/10.3389/fpls.2017.00728

BookModeling plant tissue growth and cell divisionIn: Morris R. J, editors. Mathematical Modelling in Plant Biology. Cham: Springer. pp. 107–138.https://doi.org/10.1007/9783319990705

Revisiting the Female Germline and Its Expanding ToolboxTrends in Plant Science 24:455–467.https://doi.org/10.1016/j.tplants.2019.02.003

Germline and pluripotent stem cellsCold Spring Harbor Perspectives in Biology 7:a019422.https://doi.org/10.1101/cshperspect.a019422

Discrete Laplace–Beltrami operators for shape analysis and segmentationComputers & Graphics 33:381–390.https://doi.org/10.1016/j.cag.2009.03.005

ConferenceA discussion on mixed (Longest‐Side midpoint insertion) Delaunay techniques for the triangulation refinement problemProceedings 4th International Meshing Roundtable. pp. 335–346.

Stressdependent finite growth in soft elastic tissuesJournal of Biomechanics 27:455–467.https://doi.org/10.1016/00219290(94)900213

Canalization: Genetic and Developmental AspectsAnnual Review of Ecology and Systematics 22:65–93.https://doi.org/10.1146/annurev.es.22.110191.000433

Growth and Development of ThreeDimensional Plant FormCurrent Biology 27:R910–R918.https://doi.org/10.1016/j.cub.2017.05.079

Dynamics of plant DNA replication based on PCNA visualizationScientific Reports 6:29657.https://doi.org/10.1038/srep29657
Article and author information
Author details
Funding
Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (310030L_170167)
 Célia Baroux
Agence Nationale de la Recherche (16CE930002)
 Daphné Autran
Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (310030B_160336)
 Ueli Grossniklaus
Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (IZCOZ0_182949)
 Célia Baroux
Kommission für Technologie und Innovation (16997)
 Célia Baroux
Baugarten Stiftung
 Célia Baroux
Consejo Nacional de Ciencia y Tecnología (438277)
 Elvira HernandezLagana
Forschungskredit Universitaet Zuerich (FK745020401)
 Gabriella Mosca
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank S Strauss (Max Planck Institute for Plant Breeding Research, Cologne, Germany) and B Lane (John Innes Centre, Norwich, UK) for providing code and assessment for mesh generation and simulations; S Guyer (Bitplane, Switzerland) and M Lartaud (CIRAD, France) for their help with IMARIS analyses; O Leblanc (IRD Montpellier, France) for his help with data visualization in R; the Montpellier imaging facility MRI, member of the national infrastructure FranceBioImaging supported by the French National Research Agency (ANR10INBS04) and the Center for Microscopy and Image Analysis of the University of Zurich (ZMB) for microscopy imaging infrastructures; M Ueda (Nagoya University, Japan), I Siddiqi (Center for Cellular and Molecular Biology, India), N de Storme (Katholieke Universiteit Leuven, Belgium), S Matsunaga (Tokyo Science University, Japan), M Tucker (University of Adelaide, Australia), R Dixit (Washington University in St. Louis, USA), DW Ehrhardt (Carnegie Institution for Science, Stanford, USA), WC Yang (Chinese Academy of Sciences University) and the NASC stock center for seeds. This work was funded by the University of Zürich, the IRD, a PRCI grant from the Agence Nationale de la Recherche/Swiss National Science Foundation (#ANR16CE930002/SNF310030L_170167) to DA and CB; grants from the Swiss National Science Foundation (# 310030B_160336 to UG, # IZCOZ0_182949 to CB), from the Commission for Technology and Innovation (CTI grant #16997) to CB, from the Baugarten Stiftung Zürich to CB, a CONACYT fellowship (# 438277) to EHL, and a fellowship from the Forschungskredit of the University of Zurich (#FK745020401) to GM.
Version history
 Received: December 22, 2020
 Accepted: May 3, 2021
 Accepted Manuscript published: May 7, 2021 (version 1)
 Version of Record published: June 22, 2021 (version 2)
Copyright
© 2021, HernandezLagana et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,544
 Page views

 460
 Downloads

 15
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Developmental Biology
 Immunology and Inflammation
During embryogenesis, the fetal liver becomes the main hematopoietic organ, where stem and progenitor cells as well as immature and mature immune cells form an intricate cellular network. Hematopoietic stem cells (HSCs) reside in a specialized niche, which is essential for their proliferation and differentiation. However, the cellular and molecular determinants contributing to this fetal HSC niche remain largely unknown. Macrophages are the first differentiated hematopoietic cells found in the developing liver, where they are important for fetal erythropoiesis by promoting erythrocyte maturation and phagocytosing expelled nuclei. Yet, whether macrophages play a role in fetal hematopoiesis beyond serving as a niche for maturing erythroblasts remains elusive. Here, we investigate the heterogeneity of macrophage populations in the murine fetal liver to define their specific roles during hematopoiesis. Using a singlecell omics approach combined with spatial proteomics and genetic fatemapping models, we found that fetal liver macrophages cluster into distinct yolk sacderived subpopulations and that longterm HSCs are interacting preferentially with one of the macrophage subpopulations. Fetal livers lacking macrophages show a delay in erythropoiesis and have an increased number of granulocytes, which can be attributed to transcriptional reprogramming and altered differentiation potential of longterm HSCs. Together, our data provide a detailed map of fetal liver macrophage subpopulations and implicate macrophages as part of the fetal HSC niche.

 Developmental Biology
 Neuroscience
Autism spectrum disorder (ASD) is defined by common behavioral characteristics, raising the possibility of shared pathogenic mechanisms. Yet, vast clinical and etiological heterogeneity suggests personalized phenotypes. Surprisingly, our iPSC studies find that six individuals from two distinct ASD subtypes, idiopathic and 16p11.2 deletion, have common reductions in neural precursor cell (NPC) neurite outgrowth and migration even though whole genome sequencing demonstrates no genetic overlap between the datasets. To identify signaling differences that may contribute to these developmental defects, an unbiased phospho(p)proteome screen was performed. Surprisingly despite the genetic heterogeneity, hundreds of shared ppeptides were identified between autism subtypes including the mTOR pathway. mTOR signaling alterations were confirmed in all NPCs across both ASD subtypes, and mTOR modulation rescued ASD phenotypes and reproduced autism NPCassociated phenotypes in control NPCs. Thus, our studies demonstrate that genetically distinct ASD subtypes have common defects in neurite outgrowth and migration which are driven by the shared pathogenic mechanism of mTOR signaling dysregulation.