In the SAM, DII-VENUS fluorescence reports auxin concentration with cellular resolution (Vernoux et al., 2011; Brunoud et al., 2012). To extract quantitative information about auxin distribution, we generated a DII-VENUS ratiometric variant, hereafter named qDII (quantitative DII-VENUS). qDII differs from previously used tools (Liao et al., 2015) in producing DII-VENUS and a non-degradable TagBFP reference stoichiometrically from a single RPS5A promoter (Wend et al., 2013; Goedhart et al., 2011; Figure 1—figure supplement 1A–H). By introducing a stem cell-specific pCLV3:mCherry nuclear transcriptional reporter into plants expressing qDII (Pfeiffer et al., 2016) we generated a functional and robust geometrical reference for the SAM center (Figure 1B,C and Figure 1—figure supplement 1I–M).

All analyzed meristems (21 individual SAM) showed qDII patterns similar to those obtained with DII-VENUS, with locations of auxin maxima following the phyllotactic pattern (Vernoux et al., 2011; Figure 1B–E). Despite the fact that SAMs were imaged independently and not synchronized, qDII patterns appeared highly stereotypical with easily identifiable fluorescence maxima and minima. This was confirmed by image alignment using SAM rotations (applying prior mirror symmetry if necessary; Figure 1D and Figure 1—figure supplement 2A–C). All images could be superimposed preserving the spatial distribution of auxin maxima and minima (Figure 1—figure supplement 2B). Our analysis shows that auxin distribution follows the same synchronous pattern across a population of SAMs, with low angular and rhythmic variability (Figure 1—figure supplement 2D–E, Appendix 2), with apparent stationarity up to a 137° rotation (Figure 1H).

To further quantify auxin distribution, we developed a mostly automated computational pipeline to measure SAM fluorescence (Appendix 3) (Cerutti et al., 2020). We used the spatial distribution of 1-DII-VENUS/TagBFP as a proxy for auxin distribution, hereafter named ‘auxin’ (Figure 1C) and focused on the epidermal cell layer (L1) where organ initiation takes place (Jonsson et al., 2006; Kierzkowski et al., 2013; Smith et al., 2006b; Reinhardt et al., 2003b). The location of the absolute auxin maximum value was defined as Primordium 0 (P₀). Other local maxima with lower auxin values were called P_n (Appendix 1), with n corresponding to their rank in the phyllotactic spiral (Figure 1C and Figure 1—figure supplement 2B). Note that the dynamic range of qDII allows measuring an auxin value for the vast majority of cells in the PZ and only a few cells at P0 had undetectable values of DII-VENUS, leading to an auxin value of 1. The pipeline then permits the quantification of nuclear signals and aligns all the SAMs onto a common clockwise reference frame with standardized x,y,z-orientation and with the P₀ maximum to the right. This automatic registration confirmed that auxin maxima follow a phyllotactic pattern with a divergence angle close to 137.5° (Figure 1—figure supplement 2F). It also demonstrated that maxima are positioned with a precision close to the size of a cell both in distance from the SAM center and in azimuth (angular distance) with a maximal standard deviation of 8.4 µm or 1.5 cell diameters (Figure 1G).

We then considered the temporal changes in auxin distribution by using time-lapse images over one plastochron, which corresponds to the period of this rhythmic system. P₀ and successive auxin maxima moved radially (Figure 1—figure supplement 2D). Remarkably, while the average radial distance from each local maximum P_n to the SAM center progresses (Figure 1—figure supplement 2G), the spatial deviation of this distance does not change significantly over time, reflecting the synchronized movement of local maxima, with limited meristem to meristem variation. After 10 hr, every P_n local maximum has almost reached the starting position of the next local maximum, P_n+1, but after 14 hr they have passed this position (Figure 1—figure supplement 2G). This suggests that a rotation of 137.5°, which replaces P_n by P_n+1, corresponds to a temporal progression of 10 to 14 hr (Figure 1H). This was supported by dissimilarity measurements obtained using different rotation angles between maps (Figure 1—figure supplement 2H), allowing us to confirm that plastochron last 12h ± 2h. We could thus derive a continuum of primordium development by placing P_n+1 time series one plastochron (12 hr) after P_n time series on a common developmental time axis (Figure 1H). Together with the observed developmental stationarity, this permitted the reconstruction of auxin dynamics over several plastochrons from observations spanning only one. The resulting quantitative temporal map of auxin distribution in the SAM reveals the dynamic genesis of auxin maxima in the PZ first as finger-like protrusions (visible at P_-2, P_-1 and P₀) from a permanent high auxin zone at the center of the SAM (Figure 1—figure supplement 2I–L and Video 1), as previously predicted (de Reuille et al., 2006). At later stages, auxin maxima become confined to fewer cells while auxin minima are progressively established precisely in between auxin maxima and the CZ (Figure 1—figure supplement 3).

Video 1

Download asset

We next wondered whether the motion of auxin maxima and minima could result purely from cellular growth, an hypothesis used in several theoretical models (Douady and Couder, 1996; Jonsson et al., 2006; Smith et al., 2006b; Heisler and Jönsson, 2006). By following a P₁ maximum, we observed that cells within the auxin maximum zone closest to the CZ at time 0 hr gradually transfer to the depletion zone at time 10 hr (Figure 2A–C; nuclei circled in white). At the same time, cells on the distal edge of the maximum zone show a progressive increase in their auxin level (Figure 2A–C; nuclei circled in red), suggesting a spatial shift of the auxin maximum relatively to the cellular canvas. To explore further this phenomenon, we used nuclear motion to estimate cell motion vectors and compare them with the motion of the center of auxin maximum zones, we further found that the average radial speed of auxin maxima between stages P₁ and P₄ can surpass the average displacement of individual nuclei, with a peak velocity of more than 1 µm/h at the P₂ stage (Figure 2D–E). These results show that auxin maxima are not attached to specific cells; instead they travel through the tissue, resulting in an apparent centrifugal wave of auxin accumulation. Consequently, the SAM cellular network provides a dynamic medium in which auxin maximum zones can move radially with their own apparent velocity relative to the growing tissue (Figure 2D–E). Analysis on time-courses of up to 14 hr revealed significant auxin variations in certain cells over one plastochron while auxin levels remained unchanged in others (Figure 2F–G). However, neighboring cells always showed limited differences in their temporal auxin profiles (Figure 2F–G). We concluded from these observations that there is a high definition spatio-temporal distribution of auxin, with auxin apparent movement occurring faster than growth within the tissue and providing cells with graded positional information in space and time (Figure 2H).

Figure 2

Download asset Open asset

The creation of auxin maxima first as protrusions of a high auxin zone in the CZ contrasts with the current vision of organogenesis being triggered by local auxin accumulation at the periphery of the CZ with concomitant auxin depletion around auxin maxima (Reinhardt et al., 2003a; de Reuille et al., 2006; Stoma et al., 2008; Jonsson et al., 2006; Smith et al., 2006b). This, in addition to the partial uncoupling of auxin distribution dynamics and growth, led us to reevaluate the spatio-temporal patterns of PIN1 localization, given their central role in controlling auxin distribution (Reinhardt et al., 2003a; de Reuille et al., 2006; Jonsson et al., 2006; Smith et al., 2006b). Co-visualization of a functional PIN1-GFP (Benková et al., 2003) and qDII/CLV3 fluorescence over time showed that PIN1 concentration increases from P₀ and reaches a maximum at P₂ before decreasing (Figure 3A,H and Figure 3—figure supplement 2), consistent with previous observations (Heisler et al., 2005; Bhatia et al., 2016; Caggiano et al., 2017). To quantify PIN1 cell polarities, we used confocal images after cell wall staining with the fluorescent dye propidium iodide (PI) as a reference to position the PIN1-GFP signal relative to the L1 anticlinal cell walls at each cell-cell interface (Shi et al., 2017; Figure 3B and Appendix 4). This allowed us to compute PIN1-GFP polarity for each cell-cell interface of the SAM by extracting the 3D distribution of fluorescence for PI and GFP and quantifying the difference of intensity on membranes on both sides of the cell wall (Figure 1C and Appendix 4). These cell interface polarities measure in which direction each cell interface locally contributes to orient the flow of auxin transport. Using super-resolution radial fluctuation (SRRF) microscopy (Gustafsson et al., 2016) on the same samples, we could show that this method recovers cell interface PIN1 polarities with an error below 10% (8 out of 94 interfaces analyzed). When calculating cellular PIN1 polarity vectors by integrating the cell interface polarity information for each cell, we could further show that more than 80% of the cellular polarities deviate by less than 30° between the two approaches. This quantitative evaluation (Figure 3D–G, Figure 3—figure supplement 1 and Appendix 4) validates the robustness of our method, showing that, in spite of a coarse image resolution, a vast majority of cellular polarity directions are consistent with super resolution imaging techniques. Our approach is therefore particularly suitable for monitoring global trends at the scale of a tissue. Local averaging of the cellular vectors obtained from confocal images was then used to calculate continuous PIN1 polarity vector maps in order to identify the dominant trends in auxin flux directions in the SAM (Figure 3I, and Appendix 4). At the tissue scale, the vector maps demonstrate a strong convergence of PIN1 toward the center of the SAM (Figure 3I–J and Figure 3—figure supplement 2). In addition, PIN1 polarities deviate locally toward the radial axes followed by auxin maxima when they protrude from the CZ. We detected the previously observed inversion of PIN1 polarities at organ boundaries (Heisler et al., 2005) and our quantifications show that this occurs only from P₇ (Figure 3—figure supplement 2C), thus isolating the flower from the rest of the SAM from this late stage. P₃ to P₅ show a general flux toward the SAM that is locally deflected around the zones of auxin minima before converging back toward the SAM center (Figure 3I and Figure 3—figure supplement 2). Over the course of one plastochron, only limited changes in the PIN1 polarities are observed (Figure 3—figure supplement 2), suggesting that changes in auxin distribution at this time resolution do not require major adjustments in the direction of auxin efflux at the tissue scale.

Figure 3 with 3 supplements see all

Download asset Open asset

We next asked where auxin could be produced in the SAM. YUCCAs (YUCs) have been shown to be limiting enzymes for auxin biosynthesis (Cheng et al., 2006; Liu et al., 2016). We thus mapped expression of the eleven YUC encoding genes in the SAM, using GFP reporter lines with a promoter fragment size shown to be functional for YUC1,2 and 6 (Figure 3—figure supplement 3A–N; Liu et al., 2016; Robert et al., 2013). Only YUC1,4,6 were expressed (Figure 3K, Figure 3—figure supplement 3A–F). While YUC6 showed a very weak expression in the CZ, both YUC1 and YUC4 are expressed in the L1 layer on the lateral sides of the SAM/flower boundary from P₃ for YUC4 (Figure 3K) and P₄ for YUC1 (Figure 3—figure supplement 3A and D; Cheng et al., 2006). From P₄, YUC4 expression extends over the entire epidermis of flower primordia. This is coherent with genetic and other expression data (Supplementary file 1; Cheng et al., 2006; Armezzani et al., 2018). In addition, yuc1yuc4 loss-of-function mutants show severe defects in SAM organ positioning and size (Shi et al., 2018; Pinon et al., 2013; Figure 3L and Figure 3—figure supplement 3O–U). Taken with the organization of PIN1 polarities, these results suggest that P₃-P₅ are auxin production centers for the SAM that regulate phyllotaxis and that PIN1 polarity organization allows for pumping auxin away from these production centers and towards the meristem.

In conclusion, our results suggest a scenario in which auxin distribution depends on high concentrations of auxin at the center of the SAM, and also at P_-1 and P₀, acting as flux attractors and on auxin production primarily in P₃-P₅ (Figure 3M).

To assess quantitatively whether and how the spatio-temporal distribution of auxin is interpreted in the SAM, we next introduced the synthetic auxin-induced transcriptional reporter DR5 (Friml et al., 2003; Sabatini et al., 1999; Ulmasov, 1997) driving mTurquoise2 into the qDII/CLV3 reporter line (Figure 4A–D). Cells expressing DR5 closest to the CZ were robustly positioned at an average distance of 32 µM ± 7 (SD) from the center. This corresponds to a distance at which the intensity of CLV3 reporter expression is less than 5% of its maximal value (Figure 4—figure supplement 1A). The distance from the center at which transcription can be activated by auxin is thus defined with a near-cellular precision.

Figure 4 with 1 supplement see all

Download asset Open asset

To obtain a global vision of how auxin-controlled transcription is related to auxin concentration, we performed a Principal Component Analysis (PCA) using quantified levels of DR5, auxin and CLV3 in each nucleus of the PZ during a 10 hr time series, together with their distance from the center (Figure 4E). With the first two axes accounting for around 75% of the observed variability, we unexpectedly observed orthogonality between auxin input and DR5 output, clearly marking the absence of a general correlation in the SAM (Figure 4E, inset). This unexpected finding was confirmed by the low Pearson correlation coefficients between DR5 and auxin values at the cell-level (Figure 4—figure supplement 1B). We refined our analysis by focusing on the different primordia regions. We assembled all the observed couples of values (auxin, DR5), averaged over each primordium region, on a single graph (Figure 4F). This demonstrated that, spatially, a given auxin value does not in general determine a specific DR5 value. However, values corresponding to primordia at consecutive stages follow loop-like counter-clockwise trajectories in the auxin x DR5 space (indicated by the arrow in Figure 4F). Such trajectories are symptomatic of hysteresis reflecting the dependence of a system on its history. In other words, it appears that the relationship between auxin level and DR5 expression is not direct, but is affected by another factor depending on the previous developmental trajectory of each cell (determined by parameters such as genetic activity, protein content, signal exposure, chromatin state).

We then tried to identify what in this developmental history can explain the observed differences in DR5 response to auxin. We first used our reconstructed continuum of primordium development to study the joint temporal variations of DR5 and auxin within a group of cells during primordium initiation (Appendix 5). This showed that the start of auxin-induced transcription follows the build-up of auxin concentration with a delay of nearly one plastochron (Figure 4G–H). The duration of the observed phenomenon suggests the existence of an additional process, over and above fluorescent protein maturation (Vernoux et al., 2011; Balleza et al., 2018), that creates a significant auxin response delay in primordium cells during development. Due to this delay, DR5 is not a direct readout of auxin concentration, explaining the absence of correlation between DR5 expression and auxin levels in these cells.

We next wondered what could explain a time-dependent acquisition of cell competence to respond to auxin. A first possible scenario is that cells exiting the CZ proceed through different stages of activation of an auxin-independent developmental program enabling them to sense auxin only after a temporal delay. A second possibility is that auxin controls this developmental program through a time integration process. In this scenario, cells exiting the CZ would need to be exposed to high auxin concentrations for a given time to build up an auxin transcriptional response. To test these scenarios, we treated SAMs with auxin for different periods using physiologically relevant concentrations (Reinhardt et al., 2000; Figure 5A–I). All treatments, even the shorter ones, equally degraded DII-VENUS throughout the PZ (Figure 5—figure supplement 1A–I). This suggests that auxin uptake was similar throughout the PZ, although we cannot totally discard that some differences exist. In the shorter auxin treatments (30’ and 120’), the auxin transcriptional response was mainly enhanced at P-₁ and P-₂ and to a lesser extent at the position of the predicted P-₃that is where cells are already being exposed to auxin (Figure 5I). The longer auxin treatments (300') lead to an activation of signaling in most cells in the PZ and organs, with the strongest activation being observed again at P_-1 and P_-2 but also at the predicted azimuth for P_-3, P_-4 and P_-5 (Figure 5H–I). We could further show that a 300’ treatment with a lower auxin concentration (200 nM) activated signaling similarly (at P_-1) or more strongly (at P_-2, P_-3, P_-4 and P_-5) than a 120’ 1 mM auxin treatment. Conversely, a 120’ treatment with higher auxin concentration (5 mM) lead to an activation of signaling almost as strongly as a 300’ 1 mM treatment at P_-1, although the activation was lower at P_-2 (Figure 5I). In all treatments, no significant effect was detected at P₀, consistent with the fact that DR5 activation is already maximal at this stage of development (Figure 4). We next treated pinoid (pid) mutant SAMs with exogenous auxin. pid mutants are strongly affected in polar auxin transport and in aerial organ production (Reinhardt et al., 2003a; Friml et al., 2004; Christensen et al., 2000). DR5 expression was low and radially uniform in pid SAMs, suggesting a uniform auxin distribution (Figure 5J; Friml et al., 2004). When treated with 1 mM auxin, DR5 could be activated in all cells of the periphery of the SAM (suggesting an uptake throughout the PZ as in the wild-type) only with a 300’ treatment, while a 120’ treatment had only a weak effect (Figure 5J–M and Figure 5—figure supplement 2A–C). This indicates that, even with the reduced complexity in PZ patterning of the pid mutant (Friml et al., 2004), activation of auxin signaling is still dependent on the time of exposure to auxin in all cells surrounding the CZ. Taken together, our observations support the second scenario, with the activation of signaling being a function of both time of exposure to auxin and auxin concentration. Conversely, our results are incompatible with the first scenario, where the capacity of the cells to respond to auxin is intrinsic and is not dependent upon auxin exposure time. Notably, the results with pid SAMs suggest that all cells at the SAM periphery show no intrinsic differences in their capacity to respond to auxin, in agreement with published data (Reinhardt et al., 2003a; Heisler et al., 2005; Smith et al., 2006a). Our results thus support the hypothesis that temporal integration of auxin concentration is required for downstream transcriptional activation in the SAM.

Figure 5 with 2 supplements see all

Download asset Open asset

The Auxin Response Factor (ARF) ETTIN (ETT/ARF3) plays an important role in promoting organogenesis in the SAM (Wu et al., 2015; Chung et al., 2019). Despite the fact that ETT is a non-canonical ARF, genetic data indicate that it acts together with ARF4 and MONOPTEROS/ARF5 to promote organogenesis at the SAM. We found that in a loss-of function ett3 mutant the expression of DR5 was restricted to only 2–3 cells at sites of organogenesis, an observation consistent with a role for ETT in promoting organogenesis. In addition, a 300’ 1 mM auxin treatment did not induce DR5 in the SAM (Figure 5N–Q and Figure 5—figure supplement 2D–E). Auxin signaling and ARF3 in particular have been shown to act by modifying acetylation of histones (Wu et al., 2015; Chung et al., 2019; Long et al., 2006). Pharmacological inhibition of histone deacetylases (HDACs) alone was able to trigger concomitant activation of DR5 at P₀ and P_-1 sites in the SAM (Figure 5—figure supplement 1Q–S). Taken together, these results suggest that auxin signal integration likely depends on a functional ARF-dependent auxin nuclear pathway.

Phyllotaxis is perturbed in ett mutant SAMs (Figure 5—figure supplement 1M–P; Simonini et al., 2017). Our results thus suggest that a perturbation of the temporal reading of auxin information can result in phyllotaxis defects. Supporting this idea, we also found that daily exogenous auxin treatments at the SAM affected phyllotaxis and that the efficiency of the treatment increased with both auxin concentration and treatment length. This was particularly evident for 30’ and 120’ treatments (Figure 5—figure supplement 1J–L). 300’ treatments were less efficient at higher auxin concentrations, possibly due to compensation mechanisms. These results suggest that temporal integration of auxin information at the SAM is essential for phyllotaxis.

Seeds were directly sown in soil, vernalized at 4 °C, and grown for 24 days at 21 °C under long day condition (16 hrs light, LED 150µmol/m²/s). Shoot apical meristems from inflorescence stems with a length between 0.5 and 1.5 cm were dissected and cultured in vitro as described in Prunet et al. (2016) for 16 hrs. When required, meristems were stained with 100 µM propidium iodide (PI; Merck) for 5 min. Auxin treatments were performed by immersing meristems in solutions containing indicated concentrations of indole-acetic acid (IAA) and 10 mM MES-hydrate (buffer) for indicated periods of time. Trichostatin A (TSA – Invivogen) was added to the culture medium to a final concentration of 5 µM. Meristems were cultured in TSA for 16 hrs prior to auxin treatment. For time lapses, the first image acquisition (T=0) corresponds to 2 hrs after the end of the dark period. In planta treatments were carried out on 24 day-old Col-0 plants by dropping 10 µL of IAA solution (IAA at different concentrations, 10 mM MES-hydrate and 0.01% v/v Tween-20) onto the SAM, followed by incubation for indicated lengths of time. Meristems were then washed with 100 µL of 10 mM MES buffer with 0.01% v/v Tween-20. Treatments were carried out on 5 consecutive days and perturbations in organ positioning were recorded 7 days after the end of the treatments.

Previously published transgenic lines used in this study are PIN1-GFP (Benková et al., 2003), pCLV3:mCherry-NLS (Pfeiffer et al., 2016), pYUC1-11:GFP and yuc1 yuc4/+ pDR5rev::GFP (Liu et al., 2016; Robert et al., 2013), ett-22 (Pekker et al., 2005), pid-14 (Huang et al., 2010). pRPS5a:DII-VENUS-N7-p2A-TagBFP-SV40 (qDII) and pDR5rev:2x-mTurquoise2-SV40 constructs were cloned cloned using Gateway technology (Life Sciences), and transformed in Arabidopsis thaliana (Col-0). Stable qDII homozygous lines were then crossed with pCLV3:mCherry-NLS, pDR5rev:2x-mTurquoise2-SV40 and PIN1-GFP reporter lines.

All confocal laser scanning microscopy was carried out with a Zeiss LSM 710 spectral microscope or a Zeiss LSM700 microscope. Multitrack sequential acquisitions were always performed using the same settings (PMT voltage, laser power and detection wavelengths) as follows: VENUS, excitation wavelength (ex): 514 nm, emission wavelength (em): 520–558 nm; mTurquoise2, ex: 458 nm, em: 470–510 nm; EGFP, ex: 488 nm, em: 510–558 nm; TagBFP, ex:405 nm, em: 430–460 nm; mCherry, ex: 561 nm, em: 580–640 nm; propidium iodide, ex: 488, em: 605–650 nm.

Scanning electron microscopy of meristems were carried out using a HIROX SH-3000 microscope.

Time lapses for Super Resolution Radial Fluctuation (SRRF) imaging were performed on an inverted Zeiss microscope (AxioObserver Z1, Carl Zeiss Group, http://www.zeiss.com/) equipped with a spinning disk module (CSU-W1-T3, Yokogawa, www.yokogawa.com) and a Prime95B SCMOS camera (https://www.photometrics.com) using a 63x Plan-Apochromat objective (numerical aperture 1.4, oil immersion), pixel size 175 nm or a 100x Plan-Apochromat objective (numerical aperture 1.46, oil immersion), pixel size 110 nm. GFP was excited with a 488 nm laser (150 mW) and fluorescence emission was filtered using a 525/50 nm BrightLine single-band bandpass filter (Semrock, http://www.semrock.com/). PI was excited with a 561 nm laser (80 mW) and fluorescence emission was filtered using a 609/54 nm BrightLine single-band bandpass filter (Semrock, http://www.semrock.com/). To obtain high resolution images, 200 frames were acquired with 50% laser power and 70 ms exposure time using Stream Acquisition mode. The green and red channels were acquired sequentially. For drift correction, 200 nm TetraSpeck beads (Life Technologies) were added to samples. Images were processed using the NanoJ-SRRF plugin (Gustafsson et al., 2016) with the following parameters: Ring Radius 0.5, Radiality Magnification 5, Axes in ring 6, Temporal Analysis TRPPM. SRRF time-lapses were produced by running SRRF analysis on groups of 50 frames. If aberrant PSF of Tetraspeck beads were observed, datasets were discarded.

All confocal images were pre-processed using the ImageJ software (http://rsbweb.nih.gov/ij/) for the delimitation of the region of interest. Then the CZI image files were processed using a computational pipeline relying on the numpy, scipy, pandas, czi_file Python libraries, as well as other custom libraries. Extensive details about the computational methods and algorithms are given in Appendix 3, 4 and 5.

Given the non- linear positive DR5 response, the raw values were logarithmically transformed in order to obtain a symmetric distribution of the noise. Nadaraya-Watson estimates and confidence intervals were then calculated with a confidence level of 95% in the R environment (RStudio Team, 2015). The boxplots displayed in the article were obtained by computing the median (central line), first and third quartiles (lower and upper bound of the box) and first and ninth deciles (lower and upper whiskers) using the R environment or numpy percentile function and rendered using the matplotlib Python library. Linear regressions were performed using the polyfit and polyval numpy functions. P-values were obtained using the scipy anova implementation in the f_oneway function. Principal component analysis was performed using the PCA implementation from the scikit-learn Python library. All data were generated with at least three independent sets of plants.

All experimental data and quantified data that support the findings of this study are available from the corresponding authors upon request.

Generic quantitative image and geometry analysis algorithms are provided in Python libraries timagetk, cellcomplex, tissue_nukem_3d and tissue_paredes (https://gitlab.inria.fr/mosaic/) made publicly available under the CECILL-C license. Specific SAM sequence alignment and visualization algorithms are provided in a separate project providing Python scripts to perform the complete analysis pipelines (Cerutti, 2020; copy archived at https://github.com/elifesciences-publications/sam_spaghetti).

The shoot apical meristem dynamically produces organ primordia, issuing from a central dome-shaped area, into a complex spatio-temporal pattern that is referred as phyllotaxis. In an abstract view of this structure, the meristem can be seen as dynamic collection of organ primordia characterized by their spatial trajectory relatively to the central zone (CZ) and by the evolution of their inner state. We propose a formal definition of such a system, which we name a ‘phyllotactic dynamical system’.

Let a phyllotactic dynamical system $\mathrm{S}$ be a finite set of primordia considered over a time interval $\mathrm{T}\mathrm{=}\mathrm{[}{t}_{\min }\mathrm{,}{t}_{\max }\mathrm{]}\mathrm{\subset }\mathrm{ℝ}$ and such that:

• At every time $t\in 𝒯$, each primordium $p\in 𝒮$ is characterized by its current state $\{{\tau }_p(t),x_p(t),y_p(t)\}$ where:

  • ${\tau }_p\in [{\tau }_{\min },{\tau }_{\max }]\subset \mathrm{ℝ}$ is called the developmental state of the primordium.

  • $x_p=[r_p(t),{\theta }_p(t),z_p(t)]\in {\mathrm{ℝ}}^3$ is the spatial position of the primordium in a cylindrical coordinate system, the origin of which is called the center of the system.

  • $y_p(t)\in {\mathrm{ℝ}}^d$ is a vector describing the physiological state of the primordium.

• The developmental state ${\tau }_p$ is a continuous strictly increasing function of time. Note that consequently, for every $p\in 𝒮$, ${\tau }_p$ is a bijection between $𝒯$ and ${\tau }_p\left(𝒯\right)$.

• In the case where $0\in {\tau }_p\left(𝒯\right)$, the time $t_{0,p}\stackrel{def}{=}{\tau }_p^{-1}(0)$ is called the initiation time of the primordium $p$.

• The spatial position and the physiological state are conditioned by the developmental state of the primordium in such way that:

  • (1) ${\displaystyle \{\begin{array}{ll}\mathrm{\exists }X:\left[{\tau }_{\mathrm{m}\mathrm{i}\mathrm{n}},{\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}\right]⟶{\mathbb{R}}^3& \mid \mathrm{\forall }p\in \mathcal{𝒮},\mathrm{\exists }{X}_p\in {\mathbb{R}}^3\mid \mathrm{\forall }t\in \mathcal{𝒯},x_p(t)=X_p+X\left({\tau }_p(t)\right)\\ \mathrm{\exists }Y:\left[{\tau }_{\mathrm{m}\mathrm{i}\mathrm{n}},{\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}\right]⟶{\mathbb{R}}^d& \mid \mathrm{\forall }p\in \mathcal{𝒮},\mathrm{\forall }t\in \mathcal{𝒯},y_p(t)=Y\left({\tau }_p(t)\right)\end{array}}$

• $𝒮$ is equipped with a strict total order denoted < that verifies:

  • (2) ${\displaystyle \mathrm{\forall }p,q\in \mathcal{𝒮},p<q⟹\mathrm{\forall }t\in \mathcal{𝒯},{\tau }_p(t)<{\tau }_q(t)}$

This definition reflects the idea that for any primordium, there exists an underlying physiological state, a hidden variable that determines all processes, both geometrical and physiological, that characterizes primordium development. This state can be used to rank the different organs among them, and to run through the sequence of primordia in the order of their respective development. It is actually more common to refer to primordia by their integer rank in this developmental order:

There exists a morphism between $\mathrm{(}\mathrm{S}\mathrm{,}\mathrm{<}\mathrm{)}$ and $\mathrm{(}\mathrm{\mathbb{Z}}\mathrm{,}\mathrm{<}\mathrm{)}$, and we can use it to denote the consecutiveness relationship in the strict total order of $\mathrm{S}$ by:

In a phyllotactic system, the notion of plastochron refers to the time elapsed between two consecutive organ initiations. However it is common to speak about the plastochron as a characteristic of the system when this duration does not vary over time:

We say that a phyllotactic dynamical system $\mathrm{S}$ has a plastochron $T\mathrm{\in }\mathrm{ℝ}$ if two consecutive primordia in the strict total order of $\mathrm{S}$ always have their initiation times separated by a time interval of length $T$:

A stronger assumption that is generally made on a phyllotactic system is that it develops in a steady regime, meaning that it maintains a constant rate of development. This translates into linear functions for the developmental states of primordia with a common strictly positive slope. If we add the existence of a plastochron, then this slope is naturally equal to the inverse of the plastochron:

We say that a system $\mathrm{S}$ with a plastochron $T$ has a steady development if all primordia in $\mathrm{S}$ have their developmental states increasing at the same constant rate $\mathrm{1}\mathrm{/}T$:

In such a regularly developing system, the plastochron constitutes the natural unit on the developmental scale of the primordia. Indeed, the corollary to the previous definition is that the developmental state of the primordia increases by one unit every plastochron/

In a system $\mathrm{S}$ with a steady development of plastochron $T$, all primordia increase their developmental state by 1 after a period $T$:

Another consequence is that, at the moment where a new primordium initiates, the developmental state of its immediate predecessor is exactly equal to 1. Given the steadiness of the system, this gap of one developmental unit is actually maintained throughout the evolution of primordia.

In a system $\mathrm{S}$ with a steady development of plastochron $T$, two consecutive primordia in the strict total order of $\mathrm{S}$ always have their developmental states separated by 1:

The fact that, in such a system, the primordia are all regularly staged in terms of developmental time allows to refer to them unambiguously by an integer rank from the lastly initiated primordium. Due to the previous properties, considering for each primordium the closest integer to its developmental state ensures a one-to-one mapping of primordia with a series of consecutive integers. This rank, that will remain constant during a period of one plastochron, can conversely be seen as a developmental stage through which all primordia will pass, one after the other.

In a system $\mathrm{S}$ with a steady development of plastochron $T$, at any time $t\mathrm{\in }\mathrm{T}$, the rounding function of the developmental state:

is an isomorphism. We call $\mathrm{\mathbb{P}}\mathit{}\mathrm{(}p\mathrm{,}t\mathrm{)}$ the developmental stage of primordium $p$ at time $t$. If $\mathrm{\mathbb{P}}\mathit{}\mathrm{(}p\mathrm{,}t\mathrm{)}\mathrm{=}k\mathrm{\in }\mathrm{\mathbb{Z}}$, we say that the primordium $p$ has the label ${\mathrm{\mathbb{P}}}_k$ at time t.

Intuitively, consecutive primordia should find themselves in consecutive developmental stages. The definition of developmental stage as the rounding of developmental state, and the fact that there exists a constant gap of 1 between developmental states of consecutive primordia ensures this natural property:

In a system $\mathrm{S}$ with a steady development of plastochron $T$, two consecutive primordia in the strict total order of $\mathrm{S}$ always have consecutive developmental stages:

The isomorphism between a steady phyllotactic system and a subset of consecutive integers allows to simplify the notations for the ranking of the primordia. If it was natural to denote $p+1$ the predecessor of primordium $p$, it is now possible to extend the notation to gaps of more than one unit, with an integer number that reflects the actual gap in developmental stages between two primordia.

In a system $\mathrm{S}$ with a steady development of plastochron $T$, we can extend the notation of the consecutiveness relationship in the strict total order of $\mathrm{S}$ using the isomorphism $\mathrm{\mathbb{P}}\mathit{}\mathrm{(}\mathrm{\bullet }\mathrm{,}t\mathrm{)}$ to identify the primordia by their relative developmental stages. We write this relationship as follows:

In addition to its intrinsically regular dynamics, an ideal phyllotactic system is characterized by a geometrical regularity, and the formation of spiral-like patterns. The spirals issue from the successive emergence of primordia at evenly spaced angular locations combined with an identical radial motion. In our conceptual frame, the geometrical arrangement of primordia is represented by the vectors $x_p$ of cylindrical coordinates, for which there is a common component that depends only of developmental state, and a constant part $X_p$ that depends of the considered primordium. For the system to be considered regular from a geometrical point of view, these primordium-dependent components have to follow a rigorous angular pattern, with a constant divergence angle $\alpha$ that separates two consecutive primordia.

We say that a system $\mathrm{S}$ with a steady development of plastochron $T$ has a regular phyllotaxis of divergence angle $\alpha$ if the constant parts of spatial positions of two consecutive primordia in the strict total order of $\mathrm{S}$ only differ by a rotation of angle $\alpha$ around the vertical axis:

If such a regularity property is achieved, then the system becomes highly auto-similar, so that a rotation of angle $\alpha$ corresponds to a translation in time of one plastochron $T$. In terms of primordia characteristics, this means that the spatial position $x_p$ and the physiological state $y_p$ of a primordium will be strictly identical, after one plastochron, to those of its predecessor, only up to a rotation of angle $\alpha$.

In a system $\mathrm{S}$ with a steady development of plastochron $T$ and a regular phyllotaxis of divergence angle $\alpha$, the system verifies at all times the following spatio-temporal periodicity equation:

Phyllotaxis regularity offers a way to access the ranking of primordia simply by looking at their spatial positions. If the divergence angle is such that two different primordia can not be aligned on the same direction, then the angular positions alone can be enough to provide a robust ranking of organ primordia.

In a system $\mathrm{S}$ with a steady development of plastochron $T$ and a regular phyllotaxis of divergence angle $\alpha$, if $\alpha$ is not a simple fraction of $\mathrm{2}\mathit{}\pi$, that isif:

we say that $\mathrm{S}$ has a clear regular phyllotaxis.

In a system $\mathrm{S}$ with a steady development of plastochron $T$ and a clear regular phyllotaxis of divergence angle $\alpha$, then the primordia angles $\mathrm{\{}{\theta }_p\mathrm{\mid }p\mathrm{\in }\mathrm{S}\mathrm{\}}$ are sufficient to determine the strict total order on $\mathrm{S}$:

With this conceptual frame in mind, we can define thoroughly the general problem addressed when labeling the primordia on a meristem observation, that is typically when one tries to estimate where is ${\mathrm{\mathbb{P}}}_1$ and where is ${\mathrm{\mathbb{P}}}_2$ on a microscopy image. In that problem, only a partial state is observed for each primordium, containing mostly its spatial position and possibly some quantified features. Using this information only, the goal is to stage the visible primordia, by affecting them a label that is as close as possible to their actual developmental state, in such way that if the method is used on different observations, the primordia assigned the same label really have close actual developmental states.

Given a system $\mathrm{S}$, observed at $K_t$ discrete time points $\mathrm{\{}{t}_i\mathrm{\mid }i\mathrm{\in }\mathrm{[}\mathrm{[}\mathrm{0}\mathrm{,}{K}_t\mathrm{[}\mathrm{[}\mathrm{\}}$ in which, for every $p\mathrm{\in }\mathrm{S}$ and for every $t_i$, only a partially observed state $\mathrm{\{}{\stackrel{\mathrm{~}}{x}}_p\mathit{}\mathrm{(}{t}_i\mathrm{)}\mathrm{,}{\stackrel{\mathrm{~}}{y}}_p\mathit{}\mathrm{(}{t}_i\mathrm{)}\mathrm{\}}$ is available;

Find for every $t_i\mathrm{\mid }i\mathrm{\in }\mathrm{[}\mathrm{[}\mathrm{0}\mathrm{,}{K}_t\mathrm{[}\mathrm{[}$ an estimated developmental stage function $\widehat{\mathrm{\mathbb{P}}}\mathit{}\mathrm{(}\mathrm{\bullet }\mathrm{,}{t}_i\mathrm{)}$ that verifies:

and that minimizes the average staging error of the primordia:

Now, if we make the assumption that the system shows the regularity properties detailed above (steady development and clear regular phyllotaxis of known divergence angle $\alpha$), the assignation problem can be made much simpler. Provided the system is observed over a time frame such that there is always one primordium labelled ${\mathrm{\mathbb{P}}}_k$, a solution to the problem can be found by identifying at each time point which of the primordia has its developmental stage equal to $k$. Indeed, the steady development property ensures its uniqueness, and the regular phyllotaxis allows to propagate the assignation by successive rotations of angle $\alpha$.

We say that a system $\mathrm{S}$ is ${\mathrm{\mathbb{P}}}_k$-maintaining ($k\mathrm{\in }\mathrm{\mathbb{Z}}$) if at all times, there is a primordium that has the label ${\mathrm{\mathbb{P}}}_k$:

Let $\mathrm{S}$ be a ${\mathrm{\mathbb{P}}}_k$-mantaining system with a steady development of plastochron $T$ and a clear regular phyllotaxis of divergence angle $\alpha$. The solution to the assignation of developmental stages problem can be reduced to:

Find for every $t_i\mathrm{\mid }i\mathrm{\in }\mathrm{[}\mathrm{[}\mathrm{0}\mathrm{,}{K}_t\mathrm{[}\mathrm{[}$, the primordium $p\mathrm{\in }\mathrm{S}$ such that $\widehat{\mathrm{\mathbb{P}}}\mathit{}\mathrm{(}p\mathrm{,}{t}_i\mathrm{)}\mathrm{=}k$.

The next question in order to solve this reduced problem is how to identify at each time point the primordium to label as ${\mathrm{\mathbb{P}}}_k$ based on the observations. A prerequisite to do this is that some information in the physiological state of the primordia is sufficient to know that they are currently at stage $k$. In other terms, there must be a subset of the physiological state space that is characteristic of a primordium’s developmental state $\tau$ being close to the value $k$.

Let $\mathrm{S}$ be a system with a steady development of plastochron $T$. We say that the physiological state function $Y$ is ${\mathrm{\mathbb{P}}}_k$-characteristic ($k\mathrm{\in }\mathrm{\mathbb{Z}}$) if there exists a value set ${\mathrm{\Gamma }}_k\mathrm{\subset }{\mathrm{ℝ}}^d$ such that:

In this case, the assignation of the label ${\mathrm{\mathbb{P}}}_k$ to the primordium for which the phyisiological state lies in the ${\mathrm{\mathbb{P}}}_k$-characteristic subset provides a solution to the assignation problem for which the error to minimize is, if not optimal, at least bounded by 1/2.

Let $\mathrm{S}$ be a ${\mathrm{\mathbb{P}}}_k$-mantaining system with a steady development of plastochron $T$ and a clear regular phyllotaxis of divergence angle $\alpha$ and a ${\mathrm{\mathbb{P}}}_k$-characteristic state function. The reduced solution to the assignation of developmental stages problem given by:

has a staging error ${ϵ}_{\mathrm{\mathbb{P}}}\mathit{}\mathrm{(}p\mathrm{)}$ bounded by 1/2 for every $p\mathrm{\in }\mathrm{S}$.

When the system is observed over a relatively short time period, it might be convenient to consider, for simplicity reasons, that the primordia do not change stage, and that the morphism existing between $𝒮$ and $\mathrm{\mathbb{Z}}$ at time $t_{\min }$ is preserved until $t_{\max }$. Actually, if we make this assumption for a system observed during less than one plastochron, we can show that the resulting staging error is again bounded by 1/2.

Let $\mathrm{S}$ be a system with a steady development of plastochron $T$. We say that a developmental stage assignation $\widehat{\mathrm{\mathbb{P}}}$ is stationary if it is the same at all times of observation, that is if:

If $\mathrm{S}$ is observed during a time interval $\mathrm{T}\mathrm{=}\mathrm{[}{t}_{\min }\mathrm{,}{t}_{\max }\mathrm{]}$ smaller than its plastochron $T$, then there exists a stationary developmental stage assignation with a staging error bounded by 1/2 for every $p\mathrm{\in }\mathrm{S}$:

This final observation leads us to consider that, when it applies to a system observed over less than a plastochron, the stage assignation problem has a stationary solution that guarantees an average staging error of at most 1. By labelling as ${\mathrm{\mathbb{P}}}_k$ the primordium for which the average physiological state is characteristic of stage $k$, we can obtain a staging of all primordia with the same stage at all time points without making an error of more than one developmental unit.

Let $\mathrm{S}$ be a ${\mathrm{\mathbb{P}}}_k$-mantaining system with a steady development of plastochron $T$ and a clear regular phyllotaxis of divergence angle $\alpha$ and a ${\mathrm{\mathbb{P}}}_k$-characteristic state function. If $\mathrm{S}$ is observed during a time interval $\mathrm{T}\mathrm{=}\mathrm{[}{t}_{\min }\mathrm{,}{t}_{\max }\mathrm{]}$ smaller than its plastochron $T$, then the reduced stationary solution to the assignation of developmental stages problem given by:

has a staging error ${ϵ}_{\mathrm{\mathbb{P}}}\mathit{}\mathrm{(}p\mathrm{)}$ bounded by one for every $p\mathrm{\in }\mathrm{S}$:

In a classical inhibitory field model of phyllotaxis, the developmental state ${\tau }_p=0$ of an organ primordium $p$ corresponds to the time where an initiation is decided in the peripheral zone (PZ) and would therefore match a local spatio-temporal minimum of the global inhibition field. With the idea in mind that the depletion of auxin has very often been related to the concept of ‘inhibition’ from those models of phyllotaxis, we consider that the instant where initiation happens corresponds to a local spatio-temporal maximum of auxin in the meristem. In other terms a characteristic of the primordium labelled ${\mathrm{\mathbb{P}}}_0$ should be that it has the maximal auxin level across the PZ.

Therefore if the local auxin maximality is the $j^{\mathrm{th}}$ component $Y_j\in \{0,1\}$ of the systems’s primordia state function, then we consider that the function $Y$ is ${\mathrm{\mathbb{P}}}_0$-characteristic with ${\mathrm{\Gamma }}_0=\mathrm{ℝ}\times \mathrm{\cdots }\times ]1/2,1]\times \mathrm{\cdots }\times \mathrm{ℝ}$. We observed the meristems over a time interval of 10 hr, which is less that the estimated plastochron in our experimental conditions. Therefore, by Property 9, we can define a stable assignation of developmental stages to the visible primordia of bounded error by assigning the label ${\widehat{\mathrm{\mathbb{P}}}}^{*}=0$ to the primordium that has most often the maximal value of auxin across the PZ over the times of observation. If the meristems prove to be close enough to phyllotactic systems with a plastochron and a regular phyllotaxis, then this first assignation will be enough to derive the developmental stages of all the other organ primordia based on their spatial positions. The method developed to perform this developmental stage assignation heuristic on experimental data is detailed in Appendix 2. Evidence for the regularity of the observed phyllotactic systems is discussed in Appendix 1.

In this section we develop a formal study on regularity in a phyllotactic system. Notably, we wondered to which extent the apparent similarity of the observed SAMs could be informative on the level of precision in the process of organogenesis. To answer this, we simulated a sample of phyllotactic patterns assuming that i) they are all aligned with respect to the position of their ${\mathrm{\mathbb{P}}}_0$ ii) their plastochrons and divergence angles are drawn from random distributions centered on a common average value. By varying the levels of noise on both angular positions and plastochrons, we assess how variability impacts the overlapping of phyllotactic patterns at the scale of a population.

Let us consider a 2D phyllotactic dynamical system $𝒮$ (see Appendix 1) formed by consecutive organ primordia observed on a temporal interval $𝒯$. At every time $t\in 𝒯$, each primordium labelled ${\mathrm{\mathbb{P}}}_p,p\in [[0,p_{\max }]]\subset \mathrm{\mathbb{Z}}$ is represented by its developmental stage ${\tau }_p$ and by its coordinates in a 2D cylindrical reference frame:

If we assume that the system has a plastochron $T$ and has a steady development, all primordia develop at the same constant rate $1/T$. In that case, we can derive:

In addition, if we consider that the system has a regular phyllotaxis with a divergence angle $\alpha$, that prmordia emerge on the contour of a central zone of radius $R$, and then move radially following an exponential motion law of coefficient $\beta$ we can write the state equations of the system as follows:

This can be translated into incremental equations to obtain a recursive definition of the state of the system, knowing the state of the primordium labelled ${\mathrm{\mathbb{P}}}_0$ at each time $t\in 𝒯$:

We set ourselves in a context where all the considered phyllotactic systems have previously been aligned on ${\mathrm{\mathbb{P}}}_0$, that is where $\forall t\in 𝒯,{\theta }_0(t)=0$ and ${\tau }_0(t)\in [-0.5,0.5]$.

We study what happens if we introduce variability into this system, by adding noise on two of the key variables of the system:

• A Gaussian noise of standard deviation ${\sigma }_{\alpha }$ on the divergence angle $\alpha$

• A Gaussian noise of normalized standard deviation ${\sigma }_{\tau }$ on the plastochron time $T$

To be more precise, we consider that the system still has a plastochron and a constant development rate $1/T$, but that the gap between the initiation times $t_{0,p-1}$ and $t_{0,p}$ of two consecutive primordia that should always be equal to $T$ is actually a random variable:

which translates into:

Consequently, we can formulate the recursive definition of the system as the drawing of $2p_{\max }$ random variables:

We simulate a population of such systems by generating $K_{𝒮}$ single-time instances that are all aligned on ${\mathrm{\mathbb{P}}}_0$. To do so, we draw for each instance a value for ${\tau }_0$ from a uniform distribution in [-0.5, 0.5], then use the initial values $(r_0,{\theta }_0)=(e^{\beta {\tau }_0},0)$ and construct the system recursively by drawing the corresponding random variables. This way, we obtain a population of organ primordia positions identified by their rank $p$ (Appendix 2—figure 1A).

In this random population, we are interested in which extent the generated phyllotactic patterns overlap. To do so, we estimate whether the points corresponding to primordia of the same rank can be grouped into separable clusters. Therefore, we consider the obtained primordia as a point cloud of 2D cylindrical coordinates labelled by a primordium rank:

To answer the separability question, we measure to which extent the identically labeled points ${𝒫}_p=\{(r_i,{\theta }_i),\mid p_i=p\}$ are separable by applying an unsupervised clustering algorithm. For this, we clustered the points using a k-means algorithm. The resulting clusters can be separated by linear boundaries in the Voronoi diagram associated with their centroids $\{(\widehat{{\theta }_p},\widehat{r_p})\mid p\in [[0,p_{\max }]]\}$. We measure the linear separability of primordia by looking if points with the same label gather inte the same Voronoi cell.

We use prior knowledge by setting the number of components of the k-means algorithm to $p_{\max }+1$ and by setting the initial centroids ${\displaystyle \left\{(\widehat{{\theta }_p^0},\widehat{r_p^0})\mid p\in [[0,p_{\mathrm{m}\mathrm{a}\mathrm{x}}]]\right\}}$ on the positions of the primordia in a model without any noise:

After convergence, the algorithm returns $p_{\max }+1$ centroid points that we use to construct the Voronoi diagram. This actually defines a predictor for the estimated primordium rank $\widehat{p_i}$ by looking inside which cell of the diagram lies a given point (Appendix 2—figure 1B).

Finally, we estimate the Voronoi separability $v$ of our cloud of primordia points by computing the accuracy of the primordium rank prediction (noting ${\delta }_{ij}$ the Kronecker delta on integers):

If $v$ equals to 1, it means that the primordium points group into perfectly identifiable clusters. This can be interpreted as the fact that $K_{𝒮}$ randomly sampled individuals can be superimposed perfectly, once they have been centered and aligned on their primordium which is the closest to the ${\mathrm{\mathbb{P}}}_0$ stage.

This will obviously be the case (as long as the motion coefficient $\beta$ remains realistic, typically $<1$) if no noise is introduced into the system. If ${\sigma }_{\alpha }={\sigma }_{\tau }=0$, then all individuals proceed from the same regular exact pattern, and the only variability will be the one of the instant of sampling ${\tau }_0$. We wondered up to which level of noise this separability property could be maintained, in order to understand what a high observed separability could tell us on the intrinsic regularity of a phyllotactic system.

To do so, we scanned the parameter space by varying ${\sigma }_{\alpha }$ between 0° and 20° and ${\sigma }_{\tau }$ between and 2 plastochrons, first with the values $R=30\mu m$, $\alpha ={137.51}^{\circ }$ and $\beta =0.23$ corresponding to actual measured data (Figure 1—figure supplement 2D–E). As expected, increasing the angular variability creates more elongated clusters (Appendix 2—figure 1C) that still appear separated. Yet, the Voronoi separation introduces confusion between neighboring primordia, more specifically making ${\mathrm{\mathbb{P}}}_p$ and ${\mathrm{\mathbb{P}}}_{p+3}$ overlap (Appendix 2—figure 1F). Interestingly, when we increase the plastochron variability (Appendix 2—figure 1D), the confusion concerns rather ${\mathrm{\mathbb{P}}}_p$ and ${\mathrm{\mathbb{P}}}_{p+5}$ (Appendix 2—figure 1G). In both illustrated cases, the separability score drops markedly below 95%, while the separability of the actual observed data has been evaluated in the same way at 100% (Appendix 2—figure 1E).

The landscape of separability in the ${\sigma }_{\alpha }\times {\sigma }_{\tau }$ parameter space gives an insight on the effects of variability on a population of individuals (Figure 1—figure supplement 2F–G). With no surprise, primordia points appear to be less and less identifiable as azimuthal or plastochron variability increase, and even worse when both do. But it shows that there exists a maximal level in variability up to which the clusters are still perfectly separable (Figure 1—figure supplement 2E, red contour).

The standard deviation in primordia angles measured on our observed SAMs is equal to 6.7°. We measured the resulting angular deviation in the simulated primordia points, and we could show that it depends only on the divergence angle variability ${\sigma }_{\alpha }$. Moreover, we determined the value of ${\sigma }_{\alpha }$ that matches best angular deviation of observed data for primordia ranks ranging from 1 to 5 (where the model tends to show more angular variability than the observations). This value is equal to 3.6° (Appendix 2—figure 1I).

If we fix the angular variability ${\sigma }_{\alpha }$ to this value (Figure 1—figure supplement 2E, white vertical line), then the maximal plastochron variability that is allowed for the separability to remain at 100% is close to 0.4. This means that it would be impossible to see the near-perfect superposition observed in our data if the phyllotactic system that produced it had a plastochron variability greater than 0.4, which would translate into an uncertainty on organ initiation times of nearly 5 hr. From this, we conclude that a plastochron variability of 5 hr is an upper bound of the rhythmic variability achieved by real SAMs.

Yet such a value of ${\sigma }_T$ produces primordia distributions that show much more radial variability than the observed one (Appendix 2—figure 1A vs. Appendix 2—figure 1E). To get a more precise approximation of this parameter value, we measured the resulting radial deviation in the simulated primordia points. This deviation only depends on the plastochron variability ${\sigma }_T$, and determined the value that matches best the observed radial deviations for primordia ranks ranging from 0 to 3 (where the model tends to show more radal variability than the observations). This value equals 0.22 (Appendix 2—figure 1J), which corresponds to a plausible rhythmic variability of nearly 2 hr for real SAMs. The obtained pattern is more representative of the observed deviations (Appendix 2—figure 1H) even if a more accurate model of 3D primordia distribution would be required to estimate exactly the variability parameters of the system. The approximation of 2 hr for the plastochron variability consequently validates our first order assumption that all considered SAMs are in a steady regime of development with the same plastochron duration, and that the whole set of individual primordia of a given rank forms a homogeneous population in terms of developmental state.

Another interesting feature evidenced by this analysis is the influence of the motion speed of primordia separability. We explored the $\beta \times {\sigma }_{\tau }$ parameter space by fixing the value of ${\sigma }_{\alpha }$ to the observed one, and varying the motion coefficient $\beta$ between 0 and 0.4 (Appendix 2—figure 1K). It appears that lowering the speed reduces the maximal possible value of ${\sigma }_{\alpha }$ to achieve 100% separability, as the points tend to overlap more in the radial dimension, leading to a decreasing separability at fixed angular variability.

On the other hand, increasing the speed seems to greatly affect the tolerance to plastochron variability. For instance a value of ${\sigma }_{\tau }=1$ that translates into a separability of 95% when $\beta =0.1$ suddenly drops to a separability of only 50% if $\beta$ is increased up to 0.4. However increasing motion speed does not affect this much the maximal value of ${\sigma }_{\tau }$ required to achieve 100% separability, which always remains close to 0.3. This consolidates our previous conclusions, even in the case of an underestimation of the radial speed.

Appendix 2—figure 1

Download asset Open asset

Going from microscopy images to aligned quantitative data requires a complex computational pipeline that involves several steps of image analysis, computational geometry and data manipulation. The main goal of this pipeline is to provide a representation of the signals contained in the images that allows for quantitative individual comparison and identification of invariant trends in the spatial patterns of the signals across a population. To achieve this, we need to perform three basic tasks:

• Extraction of cellular objects and signal quantification from the raw voxel intensities of the images:

  • Images are essentially structured grids with an information of signal covering a discretized space, without any explicit notion of what is an object of interest or not. It is then necessary to identify those objects within the image grid, to associate them with a spatial position and extent and to use the signal intensity information to assign quantitative values to each extracted cellular object.

• Geometrical transformation of all individual data into a common spatial reference frame:

  • In order to be able to compare individual meristems and compute statistics at the scale of a population, we need to align the spatialized data extracted from the images so that it becomes possible for instance to match organs in a comparable developmental state. To do so, we chose to use a common coordinate system into which we transform all the meristem geometries.

• Computation of a continuous representation of the signal that allows point-wise comparison:

  • The information we extract from the images is defined at the scale of cells, which provides a discrete representation of signals. If we want to compute statistics, we would need to estimate a one-to-one pairing between cells of different individuals, without being sure it exists. Instead, we decided to use a continuous representation of the signals that allows the aggregation of spatialized data. At any location in space, it would be possible to obtain a value of signal for each individual, and to compute statistics without the need for cell matching.

Appendix 3—figure 1

Download asset Open asset

The succession of computational steps that were necessary to perform these tasks is depicted in Appendix 3—figure 2. Together, they allow to reconstruct dynamic continuous representations of biological signals averaged over a population of meristems, as the maps displayed in Figure 1B, Figure 3G and I, Figure 4D and Video 1. Some tools, such as the estimation of 2D maps of epidermal signal, were used extensively throughout the analysis, even though they do not figure as boxes in the pipeline. The methods used for each of these steps are addressed in detail in the following of this Appendix.

We consider that a 3D image consists of an array $I$ of size $K_x\times K_y\times K_z$ filled with values taken in an integer intensity interval $ℐ\subset \mathrm{\mathbb{N}}$. The elements of this 3D array are called voxels. In the case of a 16-bit-encoded unsigned integer image, the intensity interval $ℐ=[[0,2^{16}-1[[$. We denote:

The voxels in the array grid can be projected into a physical space $\mathrm{\Omega }\subset {\mathrm{ℝ}}^3$ through a mapping function $\mathbf{𝐱}:[[0,K_x[[\times [[0,K_y[[\times [[0,K_z[[\to \mathrm{\Omega }$ that associates the array indices with a discrete, evenly spaced, 3D lattice ${\mathrm{\Omega }}_I\subset \mathrm{\Omega }$. The voxel size $v=(v_x,v_y,v_z)\in {\mathrm{ℝ}}^3$ defines a spacing of the lattice that is potentially different on each dimension:

This mapping allow to define the image as a function $I:{\mathrm{\Omega }}_I\to ℐ$ that associates a 3D point $\mathbf{𝐱}=(x,y,z)$ in the image definition space ${\mathrm{\Omega }}_I$ with a fluorescence intensity value $I(\mathbf{𝐱})\in ℐ$, such that $\forall (i,j,k)\in [[0,K_x[[\times [[0,K_y[[\times [[0,K_z[[$,

In the case of a multichannel image, we denote $I_S$ the image channel corresponding to the signal $S$. Nuclei are detected for each meristem acquisition independently using the pRPS5a:TagBFP (Tag) channel, which we denote $I_{\mathrm{Tag}}$.

As a first approximation, cell nuclei can be considered to appear as a roughly spherical blob of fluorescence intensity in the image $I_{\mathrm{Tag}}$, with a limitedly variable radius. To detect them we convolve the image with a sequence of $K_{\sigma }$ 3D isotropic Gaussian kernels of increasing standard deviations ${\displaystyle {\mathrm{\Omega }}_{\sigma }=\left\{{\sigma }_l\mid l\in [[0,K_{\sigma }[[\right\}}$, with ${\sigma }_0<\mathrm{\cdots }<{\sigma }_{K_{\sigma }-1}$. The response to this filtering is expected to be maximal for homogeneous spheres of intensity with a radius close to the standard deviation of the Gaussian kernel.

The sequence of filtered images can be seen as a 4D Gaussian scale-space (Lindeberg, 1994) resulting in a 4D image, which by extension we denote ${\displaystyle I_{\mathrm{T}\mathrm{a}\mathrm{g}}:{\mathrm{\Omega }}_I\times {\mathrm{\Omega }}_{\sigma }\subset {\mathbb{R}}^4\to \mathcal{ℐ}}$. In this image, the fourth dimension $\sigma$ corresponds to scale. If we note $G(\sigma )$ the discrete Gaussian kernel of standard deviation $\sigma$, the scale-space transform we use is defined by, ${\displaystyle \mathrm{\forall }\left(\mathbf{x},\sigma \right)\in {\mathrm{\Omega }}_I\times {\mathrm{\Omega }}_{\sigma }}$,

We detect nuclei as a set of local 4D response maxima in this scale-space representation. To conform with the scale-space theory, we define the scale interval ${\displaystyle {\mathrm{\Omega }}_{\sigma }}$ as a geometric sequence varying from ${\sigma }_{\min }$ to ${\sigma }_{\max }$. These two bounds are to be chosen in the typical range of variation of nuclei radius.

A 4D point $(\mathbf{𝐱},\sigma )$ is then considered a local maximum if its response $I_{\mathrm{Tag}}(\mathbf{𝐱},\sigma )$ is higher than a threshold $I_{\min }$, and higher that all the neighbouring responses at all scales. More formally, if we note $ℬ(\mathbf{𝐱},\sigma )=\{{\mathbf{𝐱}}^{\prime }\in {\mathrm{\Omega }}_I\mid \parallel {\mathbf{𝐱}}^{\prime }-\mathbf{𝐱}\parallel <\sigma \}$ the discrete ball of radius $\sigma$ in the image lattice ${\mathrm{\Omega }}_I$, the 4D point $(\mathbf{𝐱},\sigma )$ is a local maximum if it realizes the maximal value of response over ${\displaystyle \mathcal{ℬ}(\mathbf{x},\sigma )\times {\mathrm{\Omega }}_{\sigma }}$, in other terms:

Let us denote $𝒫$ the set of points $\mathbf{𝐱}\in {\mathrm{\Omega }}_I$ corresponding to a local maximum of $I_{\mathrm{Tag}}$ whose response is strictly greater than $I_{\min }$. Each point in $𝒫$ corresponds to a detected nucleus, identified by an integer index $n$ and associated with a single spatial position $P_n$ in physical coordinates, so that we can write

This detection method has been evaluated on a set of 4 manually expertized SAM images acquired at different voxel sizes with a 16 bit encoding, representing more than 5000 cells. Since the acquisition parameters for those images matched the one use in the rest of our analysis, we adjusted the method parameters to perform best on this expertized dataset. The parameter testing led to the determination of the optimal values $K_{\sigma }=3$, ${\sigma }_{\min }=0.8\mu m$, ${\sigma }_{\max }=1.4\mu m$, $I_{\min }=3000$, corresponding to an evaluated performance of 95.6% recall (percentage of the manually labelled cells that were indeed detected) and 98.5% precision (percentage of the detected cells that were actually labelled by experts).

Each image channel is quantified at the level of every detected nucleus. The signal intensity value is obtained by computing a weighted average of the channel intensity $I_S$ around the position of the nucleus. The signal images showing some local subcellular noise, the raw voxel value $I_S(P_n)$ might not be fully representative of the whole nucleus. We chose to use a distance-based Gaussian weight, of constant radius ${\sigma }_{𝒩}$ for all channels, to account for as much as possible of the signal information inside the nuclei, for which the typical measured diameter is $\sim 5\mu m$. Practically for the signal $S$, we first filter the image channel $I_S$ by a Gaussian kernel of radius ${\sigma }_{𝒩}$ and retrieve the values at the voxel positions of all detected nuclei, so that $\forall n\in 𝒩$,

For example, the local level of expression of the $CLV3$ gene, imaged using pCLV3:mCHERRY in the channel $I_{\mathrm{CLV3}}$ would we quantified as ${\mathrm{CLV3}}_n=\left(I_{\mathrm{CLV3}}*G({\sigma }_{𝒩})\right)(P_n)$. In the case of the ratiometric auxin sensor qDII, we combine the information of two fluorescence channels $I_{\mathrm{Tag}}$ and $I_{\mathrm{DII}}$ to compute the ratio of estimated signals for each nuclei point, so that $\forall n\in 𝒩$,

For the purpose of the analysis, we want to discriminate between the first layer of cells ($L_1$) and the rest of the tissue. We use an automatic method to do so, which in our case cannot rely on adjacency to the background as one would do on a segmented membrane-marker image. Instead we will use the distance of the nuclei to the estimated surface of the tissue.

This surface is computed based on the pRPSa:TagBFP channel as a 3D triangle mesh $ℳ=\{𝒱,𝒯\}$, where:

• $𝒱$ is a set of vertices

• Each vertex $v\in 𝒱$ is associated with a 3D position $M_v\in {\mathrm{ℝ}}^3$

• $𝒯$ a set of triangles defined by triplets of indices $(v_1,v_2,v_3)\in {𝒱}^3$

• $𝒯$ is such that the resulting simplicial complex forms a 2-manifold (Agoston, 2005) (Chapter 5.3: Topological Spaces, Chapter 6.3: Simplicial Complexes).

To obtain this triangle mesh, the image $I_{\mathrm{Tag}}$ is filtered by a large Gaussian kernel to diffuse nuclei intensity between cells, and is thresholded to obtain a binary region, which is meshed by applying a Marching Cubes algorithm (Lorensen and Cline, 1987) on a resampled version of the image. This mesh undergoes a phase of triangle decimation (Garland and Heckbert, 1997a) and isotropic remeshing (Botsch and Kobbelt, 2004a) to obtain a surface composed of roughly 50000 regular faces.

At this stage, we generally obtain a surface that goes both above and below the tissue (notably due to weaker intensity in the inner tissue). The remove the lower part that does not correspond to the epidermal surface, we estimate the face normal vectors to keep only the triangles for which the normal points towards the upper side of the meristem. The largest edge-connected component of this set of triangles is kept and used as the estimated meristem surface mesh $ℳ$.

With each vertex $v$ of the surface mesh $ℳ$, we associate the index $n(v)$ of the closest point in the set $𝒫$ of nuclei points: $\forall v\in 𝒱,$

We define ${ℒ}_1$ as the subset of $𝒩$ formed by indices $n$ of nuclei points $P_n$ that are the closest point to at least one vertex of $ℳ$.

The previous steps were performed individually on each frame of the time-lapse acquisitions. In our study, we focused on sequences of observation of the same individual over its development, consisting of $K_t$ multichannel images $\{I(t_i)\mid i\in [[0,K_t[[\}$, indexed by their temporal position $t_i\in \mathrm{\mathbb{N}}$ in hours relatively to the first time of acquisition $t_0=0h$. In the remaining, we will consistently index data computed from the $i$-th acquisition $I(t_i)$ by the temporal index. For instance $𝒫(t_i)=\{P_n\mid n\in 𝒩(t_i)\}$ denotes the set of nuclei points detected in $I{(t_i)}_{\mathrm{Tag}}$.

To study consistently the dynamics at the scale of the sequence of images, we need to place the quantitative nuclei information in the same spatial reference frame. To do this, we estimate 3D rigid transformations between consecutive time frames of the sequence. This estimation is performed using a block matching algorithm (Ourselin et al., 2000) applied on the pRPS5a:TagBFP channel $\{I{(t_i)}_{\mathrm{Tag}}\mid i\in [[0,K_t[[\}$ of the consecutive images. This produces $K_t-1$ isometry matrices in homogeneous coordinates $R_{t_i\leftarrow t_{i+1}}$ that can be inverted and/or multiplied to transform any frame of the sequence into the spatial reference frame of any other.

We use the registration output to transform all the detected nuclei points into the coordinate system of the first frame of the sequence. By applying the resulting rigid transforms to the the nuclei points detected at time $t_i$, we obtain a new point cloud $𝒫{(t_i)}^0$, indexed by the same set of integers $𝒩(t_i)$, such that $\forall n\in 𝒩(t_i)$,

In a second time, we want to estimate the local deformation of the tissue, notably to get a quantitative measure of individual cell motion and to approximate local cellular growth.

We compute this new transformation as a dense vector field that maps two consecutive pRPS5a:TagBFP images, that have previously been rigidly registered into the coordinate system of $I(t_0)$. We denote $\{I{(t_i)}_{\mathrm{Tag}}^0\mid i\in [[0,K_t[[\}$ the sequence of such registered images.

The vector field transforming $I{(t_i)}_{\mathrm{Tag}}^0$ into $I{(t_{i+1})}_{\mathrm{Tag}}^0$ is estimated using the block matching framework for non-linear registration (Ourselin et al., 2000). This approach has proven to be efficient on plant tissues in the case of deformations of moderate amplitude (Fernandez et al., 2010; Michelin et al., 2016) which is the case for the meristematic tissue we are considering with a $4h$ to $5h$ time interval between two frames.

The resulting vector field is actually a vectorial image ${\mathbf{𝐔}}_{t_i\to t_{i+1}}^0:{\mathrm{\Omega }}_{I(t_0)}\to {\mathrm{ℝ}}^3$ defined over the same grid as $I(t_0)$. It contains at each voxel position a 3D vector measuring the local total deformation of the tissue to go from the current frame to the next one.

We also perform the backwards non-linear registration by computing another set of non-linear transformations between $I{(t_{i+1})}_{\mathrm{Tag}}^0$ and $I{(t_i)}_{\mathrm{Tag}}^0$ in the exact same manner and obtain the vector fields ${\mathbf{𝐔}}_{t_i\leftarrow t_{i+1}}^0$ mapping the next frame to the current one.

The cellular motion between two consecutive sequence frames taken at times $t_i$ and $t_{i+1}$ in the forward direction can be estimated using the registered nuclei points $P{(t_i)}^0=\{P_n^0\mid n\in 𝒩(t_i)\}$ and the transformation ${\mathbf{𝐔}}_{t_i\to t_{i+1}}^0$ that maps the current frame into the next one by a vector field. Each nucleus $n$ can then be assigned a local displacement vector by looking at ${\mathbf{𝐔}}_{t_i\to t_{i+1}}^0(P_n^0)$.

We deduce speed vectors ${\mathbf{𝐯}}_n^0$ measuring the local speed of cellular motion by dividing the displacement vector by the time interval: $\forall n\in {ℒ}_1(t_i)$,

To infer a signal value on any 3D point in space from the values quantified at discrete nuclei points, we used an approximation strategy. Our method makes the signal continuous in space by computing a local weighted average of signal values. The weighting function $\eta$ we use is a parametric sigmoid density function of the distance $r$ to a given point, $\eta :{\mathrm{ℝ}}^{+}\to [0,1]$, which relies on two parameters: an extent parameter $R$, and a sharpness parameter $k$ such that $\eta (0)\sim 1$, $\eta (R)=\frac{1}{2}$, and $\dot{\eta }(R)=-\frac{k}{2}$ (Appendix 3—figure 3A). This sigmoid function takes the form:

The continuous signal map $\widehat{S}$ is then defined based on a point cloud $𝒫=\{P_n\in {\mathrm{ℝ}}^3\mid n\in 𝒩\}$ and the associated signal values $\{S_n\mid n\in 𝒩\}$ as : $\forall \mathbf{𝐱}\in {\mathrm{ℝ}}^3$,

Note that, for any point $\mathbf{𝐱}$ where the total density $H(\mathbf{𝐱})={\sum }_{n\in 𝒩}\eta \left(\parallel \mathbf{𝐱}-P_n\parallel \right)$ equals , the signal map is not defined. To make the map outlines closer to their actual support, we consider that $\widehat{S}(\mathbf{𝐱})$ is defined only if $H(\mathbf{𝐱})\ge \frac{1}{2}$. This constraint is equivalent to consider that the implicit surface obtained with the point cloud $𝒫$ as generator and $\eta$ as potential forms the boundary of the definition domain of the function $\widehat{S}$.

The first layer of cells forms a continuous surface and in the case of the shoot apical meristem, the curvature of this surface is such that it generally does not create overlaps in a 2D projection made along the meristem axis. Therefore, to study processes taking place at the $L_1$, we compute maps using only the subset ${ℒ}_1$ of first-layer nuclei and their 3D point cloud projected on the plane $z=0$ (Appendix 3—figure 2C-E), which we call 2D projected maps of epidermal signal. In this case, by extension we denote $\widehat{S}(x,y)=\widehat{S}(x,y,0)$.

Appendix 3—figure 2

Download asset Open asset

To determine the best parameter values $R$ and $k$ for the function $\eta$ we ran an extensive parameter exploration and measured the error made by mapping the epidermal signal to retain the values that yield the minimal error. This is measured by the average relative error between the actual signal value $S_n$ of a nucleus and the value of the 2D map $\widehat{S}$ computed with all nuclei but the considered one at the projected position of the nucleus.

We searched for optimal values on the qDII signal, as it is our main focus in this work, using the whole set of 21 available SAM image sequences. The values we obtain, shown in Appendix 3—figure 2B, are $R^{*}=7.5\mu m$ and $k^{*}=0.55\mu m^{-1}$. From now on, we will consider that, if not explicitly mentioned otherwise, the maps we use are 2D projected maps of epidermal signal computed with the optimal parameters $R=R^{*}$ and $k=k^{*}$ for the density function $\eta$.

In order to aggregate the data quantified on several individuals and expressed in their own image reference frame, we need to perform a geometrical alignment that superimposes organs with a similar developmental state.

To do so, we chose to map a common coordinate system onto the point clouds extracted from the images by landmarking a set of key geometrical features for each meristem (Appendix 3—figure 3A-B):

• the position $\mathbf{𝐜}=(x_{\mathbf{𝐜}},y_{\mathbf{𝐜}},z_{\mathbf{𝐜}})\in {\mathrm{ℝ}}^3$ of the apex center in the central zone (CZ) of the meristematic dome

• the unitary vector $\mathbf{𝐚}\in {\mathrm{ℝ}}^3$ of the main radial symmetry axis of the meristematic dome

• the unitary radial vector $\mathbf{𝐫}\in {\mathrm{ℝ}}^3$ ( $\mathbf{𝐚}⟂\mathbf{𝐫}$) of the direction of the last initiated organ primordium (labelled ${\mathrm{\mathbb{P}}}_0$) relatively to $\mathbf{𝐜}$

• the binary orientation $o\in \{-1,1\}$ of the phyllotactic spiral (clockwise or counter-clockwise)

Together, these landmarks define a new reference frame ${ℛ}^{*}$ (Appendix 3—figure 3C) into which we will rigidly transform all the nuclei points and images.

The $CLV3$ peptide, imaged in the pCLV3:mCHERRY channel, is a marker of the central zone (CZ) of the meristem that is expressed notably on the first layer of cells. We use the quantified signal values $\{{\mathrm{CLV3}}_n\mid n\in {ℒ}_1\}$ considering nuclei points from the whole sequence (registered on the reference frame ${ℛ}^0$ of the first image of the sequence) to estimate the center position.

In the resulting 2D map $\widehat{\mathrm{CLV3}}$, the CZ appears as a wide isotropic peak of signal intensity (Appendix 3—figure 3D). We assume that the CZ is the largest area of high $CLV3$ intensity in the tissue we consider. We extract a central zone domain by thresholding the 2D map $\widehat{\mathrm{CLV3}}$ and keeping the largest connected component, from which we compute the 2D center $(x_{\mathbf{𝐜}},y_{\mathbf{𝐜}})$ and the area $A_{\mathbf{𝐜}}$.

To make the estimation more robust, we perform this CZ extraction using a range of $K_{\mathbf{𝐜}}$ threshold values that depend on the average signal intensity of the sequence $\{c_k\cdot \overline{\mathrm{CLV3}}\mid k\in [[0,K_{\mathbf{𝐜}}[[\}$. In the end, we estimate $(x_{\mathbf{𝐜}},y_{\mathbf{𝐜}})$ as the average 2D center and $A_{\mathbf{𝐜}}$ and the average area of all these central zone domains obtained with $K_{\mathbf{𝐜}}=7$ and $c_k=1.2+0.1k$.

The resulting distribution of estimated CZ radii $r_{\mathbf{𝐜}}=\sqrt{\frac{A_{\mathbf{𝐜}}}{\pi }}$ among the considered individuals showed a very peaked distribution around an average value of 28 µm (Appendix 3—figure 3E).

Finally, the third coordinate $z_{\mathbf{𝐜}}$ of the 3D center of the CZ is computed as the estimated value ${\displaystyle \widehat{z^0}}$ of the $z$ coordinate of $L_1$ nuclei points in ${ℛ}^0$, taken at the center point $(x_{\mathbf{𝐜}},y_{\mathbf{𝐜}})$ detected above.

Appendix 3—figure 3

Download asset Open asset

In the reference frame of the image, the surface of the meristematic dome around the central zone might be slightly tilted. To correct this, we estimate the vector $\mathbf{𝐚}$ that corresponds to the direction around which the surface of the meristem is radially symmetrical.

This vector is characterized by a rotation matrix $R_{\mathbf{𝐚}}$ to allows to transform the reference frame ${ℛ}^0$ of the image to a new reference frame ${ℛ}_{\mathbf{𝐚}}$ where the origin is $\mathbf{𝐜}$ and the $z$ axis corresponds to $\mathbf{𝐚}$.

We estimate the rotation matrix by minimizing the tilting of the surface of the CZ in this new reference frame. To do so, we look at the dispersion of $z$ coordinates of the nuclei expressed in this reference frame as a function of their radial distance. We minimize the error between the $z$ coordinate of the $L_1$ nuclei points and their local average $\overline{z}$. To take into account the extent of the meristematic dome, this error is weighted by a Gaussian function of the distance of points to $\mathbf{𝐜}$ of standard deviation ${\sigma }_r=20\mu m$.

Note that in principle, at this stage we should re-estimate the position $\mathbf{𝐜}$ of the center using the nuclei points expressed in the reference frame ${ℛ}_{\mathbf{𝐚}}$. In our case, we considered that we could avoid this iterative estimation since the meristems are imaged from above and we can make the approximation that the vector $\mathbf{𝐚}$ forms a small angle with the $z$ axis of the images. Consequently, the initial estimation of $\mathbf{𝐜}$ in the $xy$ plane of the image does not induce a significant error.

Within the previous reference frame ${ℛ}_{\mathbf{𝐚}}$ equipped with a cylindrical coordinate system $(r,\theta ,z)$, we now look for the azimuthal direction ${\theta }_0$ of the primordium labelled ${\mathrm{\mathbb{P}}}_0$. For this, we need to locate the local spatio-temporal maximum of average auxin concentration in the peripheral zone (PZ), as explained in Appendix id1. Therefore, we look for the minimal value of the 2D map $\widehat{\mathrm{qDII}}$ (computed using the nuclei points of the whole sequence), which corresponds to a signal projected orthogonally to $\mathbf{𝐚}$.

We limit the search to the PZ by constraining that $r\in [{\rho }_{\min }{r}_{\mathbf{𝐜}},{\rho }_{\max }{r}_{\mathbf{𝐜}}]$ to avoid artifactual detections inside the CZ or in older organs (we used ${\rho }_{\min }=0.9$ and ${\rho }_{\max }=1.6$ in the analysis). In the end we get the coordinates of an auxin maximum out of which the azimuthal coordinate ${\theta }_0$ defines the unitary vector $\mathbf{𝐫}=(\cos ({\theta }_0),\sin ({\theta }_0),0)$ in the reference frame ${ℛ}_{\mathbf{𝐚}}$ (Appendix 3—figure 3F). We denote the $R_{\mathbf{𝐫}}$ the rotation matrix that allows to transform the $x$ axis of ${ℛ}_{\mathbf{𝐚}}$ into $\mathbf{𝐫}$.

Finally, rather than estimating the orientation of the meristem from the data, we chose to rely on manual expertise to determine visually wheter the arrangement of organs arond the meristem showed a clockwise ($o=-1$) or counterclockwise ($o=1$) phyllotactic spiral. This produces a simple identity or reflection matrix $R_o$, depending on the case:

In the end, the determination of the SAM landmarks $\mathbf{𝐜}$, $\mathbf{𝐚}$, $\mathbf{𝐫}$ and $o$ allows to transform the sequence registered points ${𝒫}^0$ into the common 3D reference frame ${ℛ}^{*}$ in which we will be able to compare different individuals locally (Appendix 3—figure 3G). We denote $R_{*\leftarrow t_0}$ the corresponding rigid transform, that can be written in homogeneous coordinates as

We note ${𝒫}^{*}=\{P_n^{*}\mid n\in {ℒ}_1\}$ the positions of first-layer nuclei points in this common reference frame, defined by, $\forall n\in {ℒ}_1,P_n^{*}=R_{*\leftarrow t_0}\cdot P_n^0$. From now on, we will always consider that the nuclei points have been aligned in the common SAM reference frame ${ℛ}^{*}$, and that the signal maps are computed using the aligned positions ${𝒫}^{*}$.

We also use the resulting alignment transformation to register the original images into the common reference frame by applying the same transform $R_{*\leftarrow t_0}$ to the registered image $I{(t_i)}_S^0$. This creates a registered image $I{(t_i)}_S^{*}$ expressed over a new voxel grid ${\mathrm{\Omega }}_I^{*}$ that is centered on 0 in x and y and slightly shifted in z so that:

We used the aligned 2D epidermal maps in order to create the 2D projected views of these registered images displaying only a single layer of cells that can be seen in Appendix 3—figure 3H, Figure 2A–C, Figure 3A–B, Figure 4A–C and Figure 5A–H. More specifically, we compute the values of the 2D map $\widehat{z^{*}}$ over the $(x,y)$ coordinates of the centered grid ${\mathrm{\Omega }}_I^{*}$, and produce a 2D image $I{(t_i)}_S^{*}$ defined over the same $(x,y)$ grid as the original one but keeping the voxel value of one $z$ slice per pixel, so that $\forall (x,y,z)\in {\mathrm{\Omega }}_I^{*}$,

The produced 2D image displays the intensity levels of a single curved image slice that goes through the nuclei of first cell layer, making information appear more clearly than a simple maximal intensity projection that would also include intensity from the inner layers.

In the aligned SAM reference frame ${ℛ}^{*}$, we expect to find organs in comparable developmental stages at very close spatial locations for all individuals, under the global hypothesis of a stationary and regular development. Therefore we use some a priori knowledge on regular phyllotaxis to detect the positions of the ranked organ primordia in the meristem, namely ${\mathrm{\mathbb{P}}}_0$, ${\mathrm{\mathbb{P}}}_1$, ${\mathrm{\mathbb{P}}}_2$ and so on.

Previous works on auxin dynamics in the meristem suggest that primordia correspond to local accumulation of auxin, which would be detectable as local minima in the $\widehat{\mathrm{qDII}}$ map, but also that soon after organ initiation, an auxin depletion area is formed, creating a local maximum in $\widehat{\mathrm{qDII}}$. In that respect, our primordia detection procedure consists first in detecting extremal points (both maxima and minima) in this 2D map and in labeling them in a second time by organ primordium rank.

To locate extremal regions in 2D, we are not only interested in absolute local extremality (namely peaks and troughs of the map) but also in points that are extremal in a given direction, and we therefore look for ridges and valleys of the $\widehat{\mathrm{qDII}}$ map (Appendix 3—figure 4A). Peaks, which are maximal in any direction, will appear as convergence points of ridges, troughs as convergence points of valleys, and saddle points as convergence points of a ridge and a valley.

If we consider a direction defined by a unitary vector $\mathbf{𝐮}\in {\mathrm{ℝ}}^2$, the local minimality of the 2D field $\widehat{\mathrm{qDII}}$ along $\mathbf{𝐮}$ can be defined as a scalar field ${\mathrm{\Lambda }}_{\mathbf{𝐮}}^{-}:{\mathrm{ℝ}}^2\to \{0,1\}$ that takes the value 1 only if the sign of the derivative of $\widehat{\mathrm{qDII}}$ along $\mathbf{𝐮}$ (noted ${\overrightarrow{\nabla }}_{\mathbf{𝐮}}$) changes from negative to positive. The same goes for the local maximality of the signal along $\mathbf{𝐮}$, noted ${\mathrm{\Lambda }}_{\mathbf{𝐮}}^{+}:{\mathrm{ℝ}}^2\to \{0,1\}$, except that the sign of the derivative should go from positive to negative. We define the functions ${\mathrm{\Lambda }}_{\mathbf{𝐮}}^{-}$ and ${\mathrm{\Lambda }}_{\mathbf{𝐮}}^{+}$ by: