Carlos S Galvan-Ampudia

Guillaume Cerutti

Jonathan Legrand

Romain Azaïs

Geraldine Brunoud

Steven Moussu

Christian Wenzl

Jan U Lohmann

Christophe Godin

Teva Vernoux

Romain Azais

Géraldine Brunoud

From spatio-temporal morphogenetic gradients to rhythmic patterning at the shoot apex

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

demonstrate that regulation of histone acetylation allows cells to measure the duration of the exposition to auxin preceding organ initiation, providing a remarkable example of how both spatial and temporal morphogenetic information generates rhythmic patterning.

Main text

Specification of differentiation patterns in multicellular organisms is classically thought to be regulated by gradients of morphogenetic regulators (morphogens in animals) providing positional information to cells 1 . In plants, the hormone auxin is one of the main morphogenetic regulators 2,3 . This small molecule acts during embryonic development and post-embryonic development, where it is essential for the reiterative organogenesis characteristic of plants 4 .

Notably, plant shoots develop post-embryonically through rhythmic organ generation in the shoot apical meristem (SAM), a specialized tissue with a stem cell niche in its central zone (CZ). In Arabidopsis thaliana, as in a majority of plants, organs are initiated sequentially in the SAM peripheral zone (PZ) and at a relative angle close to 137° from the previous one, either in a clockwise or anti-clockwise spiral 5 . SAM organ patterning or phyllotaxis has been extensively analyzed theoretically [START_REF] Douady | Phyllotaxis as a Dynamical Self Organizing Process Part II: The Spontaneous Formation of a Periodicity and the Coexistence of Spiral and Whorled Patterns[END_REF][7][8] . A widely accepted model proposes that the time interval between organ initiations (the plastochrone) and the spatial position of organ initiation emerge from a combined action of isotropic inhibitory signals emitted by pre-existing organs and the SAM center 5,[START_REF] Douady | Phyllotaxis as a Dynamical Self Organizing Process Part II: The Spontaneous Formation of a Periodicity and the Coexistence of Spiral and Whorled Patterns[END_REF] . Tissue growth would then self-organize organ patterning by moving organassociated signaling centers away from the stem cells and leaving space for new ones.

Biologically, evidence suggests that auxin provides the positional information that drives phyllotaxis patterning 9,10 . Auxin, thought to be synthesized throughout the meristem [11][12][13][14] , has been proposed to be transported directionally toward incipient primordia where it activates a transcriptional response leading to organ specification 3,9,15 . A network of PIN-FORMED1 (PIN1) efflux carriers, whose polarity determines the direction of auxin fluxes, regulates auxin spatio-temporal distribution cooperatively with other carriers 9,16 . The activity of this network results in accumulation of auxin that triggers organ initiation. The PIN1 network was also proposed to create an auxin depletion at the organ periphery that specifies organ boundaries and blocks organ initiation in the organ vicinity 9,[START_REF] Vernoux | The auxin signalling network translates dynamic input into robust patterning at the shoot apex[END_REF][START_REF] De Reuille | Computer simulations reveal properties of the cell-cell signaling network at the shoot apex in Arabidopsis[END_REF][START_REF] Smith | Inhibition fields for phyllotactic pattern formation: a simulation studyThis article is one of a selection of papers published on the Special Theme of Shoot Apical Meristems[END_REF][START_REF] Stoma | Flux-based transport enhancement as a plausible unifying mechanism for auxin transport in meristem development[END_REF][START_REF] Jonsson | An auxin-driven polarized transport model for phyllotaxis[END_REF] . In addition, it has been shown that the CZ is markedly less responsive to auxin [START_REF] Vernoux | The auxin signalling network translates dynamic input into robust patterning at the shoot apex[END_REF][START_REF] De Reuille | Computer simulations reveal properties of the cell-cell signaling network at the shoot apex in Arabidopsis[END_REF] . Altogether, these regional cues restrict new organ location in the growing SAM as proposed in formal models. The genetically-encoded biosensor DII-VENUS, a synthetic protein degraded directly upon sensing of auxin, recently provided a first direct experimental qualitative visualization of spatial auxin gradients in the SAM [START_REF] Vernoux | The auxin signalling network translates dynamic input into robust patterning at the shoot apex[END_REF][START_REF] Brunoud | A novel sensor to map auxin response and distribution at high spatio-temporal resolution[END_REF] .

The SAM is rather unique in that it implicates a continuous redistribution of a morphogenetic regulator in a growing tissue with helicoidal symmetry. This suggests that auxin could carry spatio-temporal morphogenetic information in the SAM. This is reminiscent of recent findings in animals that are questioning whether morphogenetic signals carry information only in space (as originally proposed in the morphogen concept [START_REF] Wolpert | Positional information and the spatial pattern of cellular differentiation[END_REF] and suggests rather a spatiotemporal nature for positional information [START_REF] Nahmad | Spatiotemporal mechanisms of morphogen gradient interpretation[END_REF][START_REF] Sagner | Morphogen interpretation: concentration, time, competence, and signaling dynamics: Morphogen interpretation[END_REF][START_REF] Wartlick | Understanding morphogenetic growth control -lessons from flies[END_REF] . Here, we used a quantitative imaging approach to reveal that auxin indeed provides spatio-temporal morphogenetic information, analyze the mechanisms generating auxin 4D dynamics and understand how this information is processed in the SAM to generate rhythmic patterning.

Spatio-temporal auxin distribution

In the SAM, the biosensor DII-VENUS fluorescence reports for auxin concentration with cellular resolution [START_REF] Vernoux | The auxin signalling network translates dynamic input into robust patterning at the shoot apex[END_REF][START_REF] Brunoud | A novel sensor to map auxin response and distribution at high spatio-temporal resolution[END_REF] . To extract quantitative data on auxin distribution, we generated a DII-VENUS ratiometric variant, named hereafter qDII (quantitative DII-VENUS). qDII consists of a RPS5A promoter driving stoichiometric co-expression of DII-VENUS and a non-degradable TagBFP reference [START_REF] Wend | A quantitative ratiometric sensor for time-resolved analysis of auxin dynamics[END_REF][START_REF] Goedhart | Quantitative co-expression of proteins at the single cell level--application to a multimeric FRET sensor[END_REF] (Extended Data Fig. 1a-h). We also introduced in plants expressing qDII a stem cell-specific CLV3::mCherry nuclear transcriptional reporter [START_REF] Pfeiffer | Integration of light and metabolic signals for stem cell activation at the shoot apical meristem[END_REF] that provided a functional and robust geometrical reference of the SAM center (Fig. 1a-b and Extended Data Fig. 1i-m).

All analyzed meristems (21 individual SAM, 24945 nuclei) showed qDII pattern similar to DII-VENUS, with auxin maxima locations following the phyllotactic pattern [START_REF] Vernoux | The auxin signalling network translates dynamic input into robust patterning at the shoot apex[END_REF] (Fig. 1a-b).

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 registration using SAM rotations (applying prior mirror symmetry if necessary; Fig 1c-d and Extended Data Fig. 2a-c). All images could be superimposed with limited loss of information definition (Extended Data Fig. 2a-c). This shows that auxin distribution follows the same synchronous pattern at population scale, with restricted angular and rhythmicity variability (Extended Data Fig. 2f-g, Supplementary Method 1), resulting, up to a rotation, in apparent stationarity (Fig. 1g).

To further quantify auxin distribution, we developed a mostly automated computational pipeline to extract SAM quantitative fluorescence (Fig. 1b and Supplementary Method 2). We used the spatial distribution of 1-DII:VENUS/TagBFP as a proxy for auxin distribution, named hereafter "auxin" (Fig. 1b) and focused on the epidermal cell layer (L1) where organ initiation takes place [START_REF] Kierzkowski | Interaction between meristem tissue layers controls phyllotaxis[END_REF] . Primordium 0 (P0) was defined as the location of the absolute auxin maximum in the PZ and other local maxima where called Pn (Supplementary Note 1), with n corresponding to their rank in the phyllotactic spiral. The pipeline then allows quantifying nuclear signal information and aligns all the SAMs to a common clockwise reference frame with standardized x,y,z-orientation with the P0 maximum to the right. This automatic registration showed that auxin maxima are positioned in the SAM with a precision close to the size of a cell both in distance to the center and in azimuth (with a maximal standard deviation of 8.4 µm or 1.5 cells, Fig. 1f).

We then considered the temporality of auxin distribution by using time-lapses over a time span close to the system period, the plastochrone. P0 and successive auxin maxima moved radially (Extended Data Fig. 2d). Remarkably, while the average radial distance of each local maximum Pn to the SAM center progresses (Extended Data Fig. 2d), their standard deviation does not change significantly in time, reflecting the synchronized movement of local maxima, with limited meristem to meristem growth variation. After 10h, every Pn local maximum has almost reached the starting position of the next local maximum, Pn+1, but after 14h they have passed it over (Extended Data Fig. 2d). This suggests that a rotation of 137.5° (Extended Data Fig. 2e) that replaces Pn by Pn+1 corresponds to a 10 to 14 hours temporal progress (Fig. 1g). This is confirmed by error measures obtained with different rotation angles (Extended Data Fig. 2h), allowing to conclude that plastochrones have a value of 12h ± 2h. We thus derived a continuum of primordium development by placing Pn+1 time series one plastochrone (12h) after Pn time series on a common developmental time axis (Fig. 1g). Together with the developmental stationarity, this allows reconstructing auxin dynamics over several plastochrones from observations spanning only one. A 4D quantitative map of auxin distribution in the SAM demonstrates dynamic building of auxin maxima in the PZ first as finger-like protrusions (visible at P-2, P-1 and P0) from a permanent high auxin zone at the center of the meristem (Supplementary video 1), as previously predicted [START_REF] De Reuille | Computer simulations reveal properties of the cell-cell signaling network at the shoot apex in Arabidopsis[END_REF] . At later stages, auxin minima are progressively established precisely in between the auxin maxima and the CZ and not surrounding the auxin maximum (Extended Data Fig. 3).

We next wondered whether the motion of auxin maxima and minima could purely result from cellular growth as proposed in previous theoretical models [START_REF] Smith | Inhibition fields for phyllotactic pattern formation: a simulation studyThis article is one of a selection of papers published on the Special Theme of Shoot Apical Meristems[END_REF][START_REF] Jonsson | An auxin-driven polarized transport model for phyllotaxis[END_REF] . Following a P1 maximum, we observed that cells in the auxin maximum at time 0h gradually became part of the depletion zone at time 10h (Fig. 2a-c). Using nuclear motion to estimate cell motion vectors and compare them with auxin maxima motion, we further found that the average radial speed of auxin maxima between stages P1 and P4 surpasses the average displacement of individual nuclei, by up to more than 1µm/h (or nearly 2 cells in 10h) at its peak in P2 stage (Fig. 2d-e).

These results show that auxin maxima are not attached to specific cells; instead they travel as a wave in the tissue. Consequently, the SAM cellular network provides a dynamic medium in which auxin distributions move radially with their own velocity relative to the growing tissue (Fig. 2d-e). Analysis on time-courses up to 14h revealed significant auxin variations in certain cells over one plastochrone while auxin levels remained unchanged in others (Fig. 2f andg, bottom rows). However, neighboring cells always showed limited differences in their temporal auxin profiles. Altogether, we concluded from these observations that there is a high definition spatio-temporal distribution of auxin that moves faster than growth in the tissue and provides cells with graded morphogenetic information both in space and time (Fig. 2h).

The control of auxin spatio-temporal dynamics

The creation of auxin maxima first as protrusions of a high auxin zone in the CZ contrasts with the largely prevailing vision of organogenesis being triggered by local auxin accumulation at the periphery of the CZ and concomitant auxin depletion around auxin maxima 9,[START_REF] De Reuille | Computer simulations reveal properties of the cell-cell signaling network at the shoot apex in Arabidopsis[END_REF][START_REF] Smith | Inhibition fields for phyllotactic pattern formation: a simulation studyThis article is one of a selection of papers published on the Special Theme of Shoot Apical Meristems[END_REF][START_REF] Stoma | Flux-based transport enhancement as a plausible unifying mechanism for auxin transport in meristem development[END_REF][START_REF] Jonsson | An auxin-driven polarized transport model for phyllotaxis[END_REF] . The partial uncoupling of auxin distribution dynamics and growth led us to question our current knowledge of the spatio-temporal patterns of PIN1, given their central role in controlling auxin distribution 9,[START_REF] De Reuille | Computer simulations reveal properties of the cell-cell signaling network at the shoot apex in Arabidopsis[END_REF][START_REF] Smith | Inhibition fields for phyllotactic pattern formation: a simulation studyThis article is one of a selection of papers published on the Special Theme of Shoot Apical Meristems[END_REF][START_REF] Jonsson | An auxin-driven polarized transport model for phyllotaxis[END_REF] . Co-visualization of a functional PIN1-GFP 3 and qDII/CLV3 fluorescence over time showed that PIN1 concentration increases from P0 and reaches a maximum at P2 before decreasing (Fig. 3a-c), consistently with previous observations 15,[START_REF] Bhatia | Auxin Acts through MONOPTEROS to Regulate Plant Cell Polarity and Pattern Phyllotaxis[END_REF] . To quantify PIN1 cellular polarities, we used staining with the fluorescent dye propidium iodide (PI) as a reference to position the PIN1-GFP signal in 3D relative to the L1 anticlinal cell walls [START_REF] Shi | Model for the role of auxin polar transport in patterning of the leaf adaxial-abaxial axis[END_REF] (Fig. 3b and Supplementary Method 2). These 3D high-resolution reconstruction of PIN1 polarities demonstrated that the crescent-shape often thought to indicate polarities in cells 15,[START_REF] Bhatia | Auxin Acts through MONOPTEROS to Regulate Plant Cell Polarity and Pattern Phyllotaxis[END_REF] does not always correlate with polarities and can thus be sometimes misleading (Extended Data Fig. 4). We mapped the vector fields of polarities to identify the trends in auxin flux directions at the scale of the SAM (Fig. 3d, and Supplementary Method 2). While some local converging PIN1 polarities can be seen close to auxin maxima at the PZ (Fig. 3e-f), the vector fields rather show large-scale convergence of PIN1 toward the center of the SAM that meet in front lines where auxin maxima protrude from the CZ (Fig. 3d and Extended Data Fig. 5). We detected the inversion of PIN1 polarities at organ boundaries previously observed 15 and our quantifications show that this occurs only from P6, thus isolating the flower from the rest of the meristem from this late stage. P3 to P5 stages show a general flux toward the SAM that is locally deflected around the zones of auxin minima before converging back toward the meristem center (Fig. 3d and Extended Data Fig. 5). Over the course of a plastochrone, only limited changes of the PIN1 network are observed (Extended Data Fig. 5), suggesting that changes in auxin distribution at this time resolution might not require important adjustments of the direction of the auxin fluxes. Other carriers are active in the SAM 16 and the role of PIN1 in auxin distribution could also have been overestimated.

The analysis of the transport network led us to question where auxin could be produced in the SAM. YUCCAs (YUCs) have been shown to be limiting enzymes for auxin biosynthesis during development 11,[START_REF] Liu | Local Transcriptional Control of YUCCA Regulates Auxin Promoted Root-Growth Inhibition in Response to Aluminium Stress in Arabidopsis[END_REF] . We thus mapped expression of the eleven YUCs in the SAM (Extended Data Fig. 6) [START_REF] Liu | Local Transcriptional Control of YUCCA Regulates Auxin Promoted Root-Growth Inhibition in Response to Aluminium Stress in Arabidopsis[END_REF][START_REF] Robert | Local Auxin Sources Orient the Apical-Basal Axis in Arabidopsis Embryos[END_REF] . Only YUC1,4,6 were detected and only YUC6 in the SAM proper with a very weak expression in the CZ (Fig. 3g-h and Extended Data Fig. 6). This is coherent with genetic and expression data (Extended Data Table 1) 11,[START_REF] Armezzani | Transcriptional induction of cell wall remodelling genes is coupled to microtubule-driven growth isotropy at the shoot apex in Arabidopsis[END_REF] . Both YUC1 and YUC4 are expressed at the L1 layer on the lateral sides of the SAM/flower boundary from P3 for YUC4 and in P4 for YUC1 11 (Fig. 3g-h and Extended Data Fig. 6). From P4, YUC4 expression extends over the entire epidermis of flower primordium. Together with the PIN1 network organization, these expression patterns thus identify the primordia at P3-P5 stages as auxin production centers for the SAM.

In conclusion, there is a structured global organization of the PIN1 pump network, with high concentrations of auxin at the center of the SAM but also at P-1 and P0, all acting as flux attractors (Fig. 3i). The stem cell niche could thus play the role of a system-wide organizer of the auxin transport network, coherently with previous observations [START_REF] De Reuille | Computer simulations reveal properties of the cell-cell signaling network at the shoot apex in Arabidopsis[END_REF] . Auxin maxima could then first emerge from the CZ along lines where auxin fluxes converge. Convergence of the fluxes towards the auxin maxima would at the same time diverge fluxes away from areas where auxin minima appear (Fig. 3i). Our data further suggest that early flowers could provide a memory of the developmental pattern and allow self-sustained distribution of auxin by redistributing auxin back in the system through biosynthesis (Fig. 3i).

The role of time in transcriptional responses to auxin

To assess quantitatively whether and how auxin spatio-temporal distribution is interpreted in the SAM, we next introduced the synthetic auxin-induced transcriptional reporter DR5 [START_REF] Sabatini | An auxin-dependent distal organizer of pattern and polarity in the Arabidopsis root[END_REF][START_REF] Ulmasov | Aux/IAA Proteins Repress Expression of Reporter Genes Containing Natural and Highly Active Synthetic Auxin Response Elements[END_REF] driving mTurquoise into the qDII/CLV3 reporter line (Fig. 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 CLV3 reporter expression is lower than 5% of its maximal value (Extended Data Fig. 7a). The distance from the center at which transcription can be activated by auxin is thus defined with a near-cellular precision.

To obtain a global vision of how auxin 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 of 10h time series together with their distance from the center (Fig. 4e). With the first two axis reporting for around 75% of the observed variability, we unexpectedly observed orthogonality between auxin input and DR5 output clearly marking an absence of a general correlation in the meristem (Fig. 4e, inset). This unexpected finding was confirmed by the low numerical values of Pearson correlation coefficients between DR5 and auxin values at cell-level (Extended Data Fig. 7b). We refined our analysis by studying locally correlations between auxin and DR5 in the different primordia regions (Fig. 4d). We observed homogeneous auxin and DR5 behaviors that characterized each region (Extended Data Table 2), showing that the spatial position of a cell conditions the link between its auxin input and the state of its transcriptional response. We then assembled on a single graph all the observed couples of values (auxin, DR5) averaged over each primordium region (Fig. 4h). This evidenced that, spatially, a given auxin value does not in general determine a specific DR5 value. However, values corresponding to primordia at consecutive stages follow looplike counter-clockwise trajectories in the auxin x DR5 space (indicated by an arrow in Fig. 4h). Such trajectories are symptomatic of hysteresis reflecting the dependence of a system on its history. In other words, DR5 expression is dependent on the developmental history of the cells.

We then used our reconstructed continuum of primordium development to study the joint temporal variations of DR5 and auxin within a local group of cells during initiation (Supplementary Method 3). This showed that the start of auxin-induced transcription follows the building-up of auxin concentration with a delay of nearly one plastochrone (Fig. 4f-g).

The duration of the observed phenomenon suggests the existence of an additional process, different from protein maturation [START_REF] Balleza | Systematic characterization of maturation time of fluorescent proteins in living cells[END_REF] , that creates a significant DR5 response delay of primordia cells to auxin during development. Due to this delay, DR5 is not a direct readout of auxin concentration, potentially explaining the absence of correlation between DR5 expression and auxin levels.

We 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 after a temporal delay. A second possibility is that auxin controls this developmental program through a time integration process. In that case, cells exiting the CZ would need to be exposed to high auxin concentration for a certain amount of time to mount up an auxin transcriptional response. To test these scenarios, we treated SAMs with auxin for different durations (Fig. 5ai). In the shorter auxin treatments (30' and 120'), auxin output was only enhanced at P0, P1 and P2 i.e. where cells have already been exposed to auxin (Fig. 5f-g and i, Extended Data Fig. 8a). On the other hand, the longer auxin treatments (300') activated signaling in most cells in the PZ and organs (Fig. 5h,i). Both observations are compatible with the second scenario where longer exposure allows activating signaling in more cells and are incompatible with the first one, where the capacity of the cells to respond to auxin is intrinsic and does not dependent upon auxin exposure time. Our results thus indicate that temporal integration of the auxin signal controls the activation of transcription in the SAM.

Chromatin state is one mechanism that allows for temporal integration of signals [START_REF] Stoma | Flux-based transport enhancement as a plausible unifying mechanism for auxin transport in meristem development[END_REF][START_REF] Nahmad | Spatiotemporal mechanisms of morphogen gradient interpretation[END_REF][START_REF] Angel | A Polycomb-based switch underlying quantitative epigenetic memory[END_REF][START_REF] Coda | Distinct modes of SMAD2 chromatin binding and remodeling shape the transcriptional response to NODAL/Activin signaling[END_REF] .

Also, ARFs and Aux/IAAs have been shown to act by modifying acetylation of the chromatin [START_REF] Martin-Arevalillo | Structure of the Arabidopsis TOPLESS corepressor provides insight into the evolution of transcriptional repression[END_REF][START_REF] Wu | Auxin-regulated chromatin switch directs acquisition of flower primordium founder fate[END_REF] . Pharmacological inhibition of histone deacetylases (HDACs) alone was able to trigger concomitant activation of DR5 at P0 and P-1 sites in the PZ (Fig. 5j-l, Extended Data Fig. 8b). This result demonstrates that the timing of auxin signaling induction at the P-1 site depends on the chromatin status. It further suggests that chromatin acetylation status of cells at the boundary of the CZ acts as a memory of their exposition time to auxin, providing a mechanism for temporal integration of auxin-based positional information.

Discussion

Through a fully quantitative analysis of auxin dynamics, our study provides a demonstration that rhythmic organ initiation at the SAM is driven by a combination of high-precision spatiotemporal graded distributions of auxin and of the use of the duration of cell exposition to auxin to differentiate temporally sites of organ initiation (Fig. 6). Such a mechanism is likely essential for rhythmic organ patterning in the SAM as auxin-based spatial information prespecifies several sites of organ initiation and is thus insufficient (Supplementary video 1).

Temporal integration of the auxin signal could occur through a chromatin acetylation mechanism. As chromatin acetylation represses auxin signaling in the CZ (Ma, Y., Miotk, A., Sutikovic, Z., Medzihradszky, A., Wenzl, C., Ermakova, O., Gaillochet, C., Forner, J., Utan, G., Brackmann, K., Galvan-Ampudia, C. S., Vernoux, T., Thomas, G. & Lohmann, J. U.

WUSCHEL acts as a rheostat on the auxin pathway to maintain apical stem cells in Arabidopsis. bioRxiv 468421; doi: https://doi.org/10.1101/468421 (2018)), this provides an original mechanism tightly linking stem cell maintenance and differentiation by precisely positioning organ initiation at the boundary of the stem cell niche while allowing for sequential organ initiation.

The existence of high definition spatio-temporal auxin gradients suggests that similarly to several morphogens in animals [START_REF] Nahmad | Spatiotemporal mechanisms of morphogen gradient interpretation[END_REF][START_REF] Dessaud | Interpretation of the sonic hedgehog morphogen gradient by a temporal adaptation mechanism[END_REF][START_REF] Scherz | Extended exposure to Sonic hedgehog is required for patterning the posterior digits of the vertebrate limb[END_REF][START_REF] Maden | Retinoid signalling in the development of the central nervous system[END_REF] the robustness of SAM patterning is at least in part due to highly reproducible spatio-temporal signal distribution. Our analysis questions how auxin transport could generate this high definition signal distribution and points again to the stem cell niche that could be crucial for organizing the transport network. The unique quantitative resource we have generated will allow addressing these questions in depth.

Methods Summary

Plant material and growth conditions

Seeds were directly sown in soil, vernalized at 4 °C, and growth for 24 days at 21 °C under long day photoperiod (16 hrs light, LED 150µmol/m²/s). Shoot apical meristems from inflorescence stems with a length between 0.5 and 1.5 cm where dissected and cultured in vitro as described in [START_REF] Prunet | Live confocal imaging of Arabidopsis flower buds[END_REF] for 16 hrs. When required, meristems were stained with 100 µM propidium iodide (Sigma) for 5 min. Auxin treatments were performed by immersing meristems with 1 mM Indole-acetic acid (IAA) for different times. Trichostatin A (TSA -Invivogen) was added to the ACM plates to a final concentration of 5 µM. Meristems were cultured in TSA for 16 hrs prior auxin treatment. For time lapses, first image acquisition (T=0) correspond to 2 hrs after the dark period, with the exception of auxin treatments where imaging is right after lights were on.

Previously published transgenic lines used in this study are PIN1-GFP 3 , promCLV3:mCherry-NLS [START_REF] Pfeiffer | Integration of light and metabolic signals for stem cell activation at the shoot apical meristem[END_REF] , and promYUC1 to 11-GFP [START_REF] Liu | Local Transcriptional Control of YUCCA Regulates Auxin Promoted Root-Growth Inhibition in Response to Aluminium Stress in Arabidopsis[END_REF][START_REF] Robert | Local Auxin Sources Orient the Apical-Basal Axis in Arabidopsis Embryos[END_REF] . Quantitative DII (qDII), promRPS5a:DII-VENUS-N7-p2A-TagBFP-SV40, reporter line and the auxin synthetic promoter DR5rev:2x-mTurqouise2-SV40 was cloned using Gateway technology (Life Sciences), and transformed in Arabidopsis thaliana (Col-0). Stable qDII homozygous lines were then crossed with promCLV3, promDR5rev:2x-mTurqouise2-SV40 and PIN1-GFP reporter lines.

Imaging

All confocal laser scanning microscopy was done with a Zeiss LSM 710 spectral microscope.

Multitrack sequential acquisition was performed using always 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-663 nm.

Quantitative image analysis

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 in a computational pipeline developed by the authors and relying essentially on the numpy, scipy, pandas, czi_file Python libraries, as well as other custom libraries. Extensive details about the developed methods and algorithms are given in Supplementary Method 2.

Statistics and reproducibility

Confidence intervals were calculated with a confidence level of 95% in the R environment [START_REF]RStudio: Integrated Development for R[END_REF] .

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 3 independent sets of plants.

Correspondence and request for materials should be addressed to T.V. (teva.vernoux@enslyon.fr) or C.G. (christophe.godin@ens-lyon.fr).

Figure 6. Spatio-temporal gradients of auxin translate into rhythmic organ patterning through time integration

A maximum of auxin protrudes from a high auxin concentration zone at the CZ faster than the cells move radially, possibly as a result of centripetal auxin fluxes. Cells exiting the CZ that are exposed to high auxin levels progressively acquire competence for transcriptional response through chromatinmediated time-integration of the auxin signal. This leads to activation of transcriptional responses with a delay close to the system period, the plastochrone. Polarity vectors detected by our method. Arrow with a green outline represents the polarity vector based on signal quantification in (g). i. Projection of PIN-signal. The arrow indicates cell 1 polarity assigned by visual impression as in b. Scale bars = 2 µM. j-m. Z-slices of image in (i) with PIN1-GFP (green) and PI (magenta) channels showing full overlap of both signals at the anticlinal cell wall between cell 1 and 2. n. Signal distribution for both signals in the wall between cell 1 and 2 as in (g). o. Polarity vectors detected. Arrow with a green outline represents the polarity vector based on signal quantification in (i).

Extended Data Figure 5. Vector fields and PIN global polarities

AT1G04180 YUC9 0 0 0 0 AT2G33230 YUC7 0 0 0 0 AT5G43890 YUC5 0 0 0 0 AT1G48910 YUC10 0 0 0 0 AT1G21430 YUC11 0 0 0 0 Transcripts Per Kilobase Million (TPM). Data extracted from 35 .

Extended Data Table 2: Auxin / DR5 variations in a primordium show a stagecharacteristic behaviour

Values of auxin and DR5 were quantified in tissue areas moving accordingly to cellular motion, across 10h of observation (rows 1 and 2) and their temporal derivatives estimated as the slope of a linear regression (rows 3 and 4). At the first order, each developmental stage from P-3 to P5 can be identified by a unique combination of the values of these four indicators. For instance a high stable value of both auxin and DR5 is characteristic of stage P0, whereas high increasing auxin with low increasing DR5 is characteristic of stage P1. "-", "0","+" and "++" correspond to arbitrary thresholds for the different variables: Auxin: "0" < 0.6 < "+" < 0.85 < "++"; DR5: "0"<5000<"+"; Auxin variation: "-"<-0.05<"0"<0.05<"+"; DR5 variation: "-"<-100<"0"<100<"+".

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".

Definition 1: Let a phyllotactic dynamical system 𝓢 be a finite set of primordia considered on a time interval 𝓣 = [𝒕 min , 𝒕 max ] ⊂ ℝ and such that:


o ∃𝑋 ∶ [𝜏 min , 𝜏 max ] → ℝ 3 ∶ ∀𝑝 ∈ 𝒮, ∃𝑋 𝑝 ∈ ℝ 3 ∶ ∀𝑡 ∈ 𝒯, 𝑥 𝑝 (𝑡) = 𝑋 𝑝 + 𝑋 (𝜏 𝑝 (𝑡))
o ∃𝑌 ∶ [𝜏 min , 𝜏 max ] → ℝ 𝑑 ∶ ∀𝑝 ∈ 𝒮, ∀𝑡 ∈ 𝒯, 𝑦 𝑝 (𝑡) = 𝑌 (𝜏 𝑝 (𝑡))  𝒮 is equipped with a strict total order denoted < that verifies:

∀𝑝, 𝑞 ∈ 𝒮, 𝑝 < 𝑞 ⇒ ∀𝑡 ∈ 𝒯, 𝜏 𝑝 (𝑡) < 𝜏 𝑞 (𝑡)

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 characterize 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:

Property 1: There exists a morphism between (𝒮, <) and (ℤ, <), and we can use it to denote the consecutiveness relationship in the strict total order of 𝒮 by:

∄𝑟 ∈ 𝒮 ∶ 𝑝 < 𝑟 < 𝑞 def ⇔ 𝑞 = 𝑝 + 1
In a phyllotactic system, the notion of plastochrone refers to the time elapsed between two consecutive organ initiations. However it is common to speak about the plastochrone as a characteristic of the system when this duration does not vary over time:

Definition 2:
We say that a system 𝒮 has a plastochrone 𝑇 ∈ ℝ if two consecutive primordia in the strict total order of 𝒮 always have their initiation times separated by a time interval of length 𝑇:

∀𝑝, 𝑝 + 1 ∈ 𝒮, 𝑡 0 𝑝+1 = 𝑡 0 𝑝 -𝑇
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:

Definition 3:
We say that a system 𝒮 with a plastochrone 𝑇 has a steady development if all primordia in 𝒮 have their developmental states increasing at the same constant rate 1/𝑇:

∀𝑝 ∈ 𝒮, ∀𝑡 ∈ 𝒯, 𝜏 𝑝 (𝑡) = 1 𝑇 (𝑡 -𝑡 0 𝑝 )
Property 2: In a system 𝒮 with a steady development of plastochrone 𝑇, all primordia increase their developmental state by 1 after a period 𝑇:

∀𝑝 ∈ ∀𝑡 ∈ 𝒯, 𝜏 𝑝 (𝑡 + 𝑇) = 𝜏 𝑝 (𝑡) + 1
Property 3: In a system 𝒮 with a steady development of plastochrone 𝑇, two consecutive primordia in the strict total order of 𝒮 always have their developmental states separated by 1: Phyllotaxis regularity offers a way to access the ranking of primordia simply by looking at their spatial positions. For instance, if the distance of a primordium to the center is increasing with time, then the order on distances reflects the order on primordia developmental states. Cues from distances and angular positions can typically be combined for an even more robust ranking of organ primordia. With this conceptual frame in mind, we can define thoroughly the problem addressed when labeling the primordia on a meristem observation, i.e. typically when one tries to estimate where is ℙ 1 and where is ℙ 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. and that minimizes the average staging error of the primordia:

∀𝑝, 𝑝 + 1 ∈ 𝒮,
𝜀 𝓈 = 1 |𝒮| ∑ 𝜀 𝓈 (𝑝) 𝑝∈𝒮 = 1 |𝒮| ∑ ( 1 𝑛 ∑|𝓈̂𝒕 𝒊 (𝑝) -𝜏 𝑝 (𝑡 𝑖 ) | 𝑛 𝑖=1 ) 𝑝∈𝒮 Definition 5:
We say that a system 𝒮 is ℙ 𝑘 -maintaining (𝑘 ∈ ℤ) if at all times, there is a primordium that has the label ℙ 𝑘 :

∀𝑡 ∈ 𝒯, ∃𝑝 ∈ 𝒮, 𝓈 𝑡 (𝑝) = 𝑘 Property 9: Reduced ℙ 𝒌 assignation problem in a clear regular phyllotaxis: Let 𝒮 be a ℙ 𝑘mantaining system with a steady development of plastochrone 𝑇 and a clear regular phyllotaxis of divergence angle 𝛼. The solution to the assignation of developmental stages problem can be reduced to:

Find for every 𝑡 𝑖 ∈ {𝑡 1 , … , 𝑡 𝑛 }, the primordium 𝑝 ∈ 𝒮 such that 𝓈̂𝒕 𝒊 (𝑝) = 𝑘.

Definition 6:

Let 𝒮 be a system with a steady development of plastochrone 𝑇. We say that the physiological state function 𝑌 is ℙ 𝒌 -characteristic (𝑘 ∈ ℤ) if there exists a value set Γ 𝑘 ⊂ ℝ 𝑑 such that:

∀𝜏 ∈ [𝜏 min , 𝜏 max ], 𝑌(𝜏) ∈ Γ 𝑘 ⇔ |𝜏 -𝑘 | < 1 2
Property 10: Developmental stage stability condition: Let 𝒮 be a system with a steady development of plastochrone 𝑇. We say that a developmental stage assignation is stable if it is the same at all times of observation, i.e. if:

∃𝓈 ∶ 𝒮 → ℤ, ∀𝑡 ∈ 𝒯, ∀𝑝 ∈ 𝒮, 𝓈 𝑡 (𝑝) = 𝓈(𝑝)
If 𝒮 is observed during a time interval smaller than its plastochron, then there exists a stable developmental stage assignation with a staging error bounded by 1/2 for every 𝑝 ∈ 𝒮:

𝑡 max -𝑡 min < 𝑇 ⇒ ∃𝓈 ∶ 𝒮 → ℤ, ∀𝑝 ∈ 𝒮, 𝜀 𝓈 (𝑝) = 1 𝑡 max -𝑡 min ∫ |𝓈(𝑝) -𝜏 𝑝 (𝑡) |𝑑𝑡 𝑡 max 𝑡 min < 1 2
Property 11: Stable ℙ 𝒌 -characterization solution: Let 𝒮 be a ℙ 𝑘 -mantaining system with a steady development of plastochrone 𝑇, a clear regular phyllotaxis of divergence angle 𝛼 and a ℙ 𝑘characteristic state function. If 𝒮 is observed during a time interval smaller than its plastochrone 𝑇, then the reduced stable solution to the assignation of developmental stages problem given by:

𝓈(𝑝) = 𝑘 ⇔ 1 𝑛 ∑ 𝑦 ̃𝑝(𝑡 𝑖 ) 𝑛 𝑖=1 ∈ Γ 𝑘
has a staging error bounded by 1 for every 𝑝 ∈ 𝒮:

𝑡 n -𝑡 1 < 𝑇 ⇒ ∀𝑝 ∈ 𝒮, 𝜀 𝓈 (𝑝) < 1
In a classical inhibitory field model of phyllotaxis, the developmental state 𝜏 𝑝 = 0 of an organ primordium 𝑝 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 0 primordium should be that it has the maximal auxin level across the PZ.

Therefore if the local auxin maximality is the 𝑗 th component 𝑌 𝑗 ∈ {0,1} of the state function, then we consider that the function 𝑌 is ℙ 0 -characteristic with Γ 0 = ℝ × … ×]1/2 , 1] × … × ℝ. We observed the meristems over a time interval of 10 hours, which is less that the estimated plastochrone in our experimental conditions. by Property 11, we can define a stable assignation of developmental stages to the visible primordia of bounded error by assigning the label 𝓈(𝑝) = 0 to the primordium that has the most often the maximal value of auxin across the PZ over the time of observation. If the meristems prove to be close enough to phyllotactic systems with a plastochrone 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 Supplementary Method 2. Evidence for the regularity of the observed phyllotactic systems is discussed in Supplementary Method 1.

Supplementary Method 1: Effects of variability on a computational phyllotactic system

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 ℙ 0 ii) their plastochrones 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 plastochrones, 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 Supplementary Note 1) formed by 𝑁 organ primordia observed on a temporal interval 𝒯. At every time 𝑡 ∈ 𝒯, each primordium labeled 𝑝 is represented by its developmental age 𝜏 𝑝 and by its coordinates in a 2D cylindrical reference frame:

𝒮 = {ℙ 𝑝 (𝑡) = (𝜏 𝑝 (𝑡), 𝑟 𝑝 (𝑡), 𝜃 𝑝 (𝑡)) ∶ 𝑝 ∈ ⟦0, 𝑁⟦}
If we assume that the system has a plastochrone and has a steady development, all primordia develop at the same constant rate 1/𝑇. In that case, we can derive:

∀𝑡 ∈ 𝒯, ∀𝑝 ∈ ⟦1, 𝑁⟦, 𝜏 𝑝 (𝑡) -𝜏 𝑝-1 (𝑡) = 1
In addition, if we consider that the system has a regular phyllotaxis with a divergence angle 𝛼, a central zone radius 𝑅, and an exponential radial motion law of coefficient 𝛽 we can write the state of the system as follows:

∀𝑡 ∈ 𝒯, ∀𝑝 ∈ ⟦0, 𝑁⟦, { 𝜃 𝑝 (𝑡) = 𝜃 0 + 𝑝 ⋅ 𝛼 𝑟 𝑝 (𝑡) = 𝑅 ⋅ 𝑒 𝛽𝜏 𝑝 (𝑡)
This can be translated into incremental equations to obtain a recursive definition of the state of the system, knowing ℙ 0 (𝑡):

∀𝑡 ∈ 𝒯, ∀𝑝 ∈ ⟦1, 𝑁⟦, { 𝜃 𝑝 (𝑡) = 𝜃 𝑝-1 (𝑡) + 𝛼 𝑟 𝑝 (𝑡) = 𝑟 𝑝-1 (𝑡)𝑒 𝛽
We set ourselves in a context where all the considered phyllotactic systems have previously been aligned on ℙ 0 , i.e. where ∀𝑡 ∈ 𝒯, 𝜃 0 (𝑡) = 0 and 𝜏 𝑝 (𝑡) ∈ [-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 𝜎 𝛼 on the divergence angle 𝛼  A Gaussian noise of normalized standard deviation 𝜎 𝑇 on the plastochrone time 𝑇

To be more precise, we consider that the system still has a plastochrone and a constant development rate 1/𝑇, but that the gap between the initiation times 𝑡 0 𝑝 and 𝑡 0 𝑝-1 of two consecutive primordia that should always be equal to 𝑇 is actually a random variable:

𝑡 0 𝑝-1 -𝑡 0 𝑝 ∼ 𝑇 ⋅ 𝒩(1, 𝜎 𝑇 )
which translates into:

∀𝑡 ∈ 𝒯, 𝜏 𝑝 (𝑡) -𝜏 𝑝-1 (𝑡) ∼ 𝒩(1, 𝜎 𝑇 )
Consequently, we can formulate the recursive definition of the system as the drawing of 2(𝑁 -1) random variables:

∀𝑡 ∈ 𝒯, ∀𝑝 ∈ ⟦1, 𝑁⟦, { 𝜃 𝑝 (𝑡) -𝜃 𝑝-1 (𝑡) ∼ 𝒩(𝛼, 𝜎 𝛼 ) 𝑟 𝑝 (𝑡)/𝑟 𝑝-1 (𝑡) ∼ 𝑒 𝛽⋅𝒩(1,𝜎 𝑇 )
We simulate a population of such systems by generating 𝐾 single-time instances that are all aligned on ℙ 0 . To do so, we draw for each instance a value for 𝜏 0 from a uniform distribution in [-0.5, 0.5], then use the initial values (𝑟 0 , 𝜃 0 ) = (𝑒 𝛽𝜏 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 𝑝 (Supplementary 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 labeled by primordium rank:

{((𝑟 𝑖 , 𝜃 𝑖 ), 𝑝 𝑖 ) ∶ 𝑖 ∈ ⟦1, 𝑁 * 𝐾⟧}
To answer the separability question, we measure to which extent the identically labeled points {(𝑟 𝑖 , 𝜃 𝑖 ) ∶ 𝑝 𝑖 = 𝑝} 𝑝∈⟦0,𝑁⟦ are separable by applying an unsupervised clustering algorithm. We chose to use the k-means method because it outputs a Voronoi diagram of the obtained centroid points, which places linear boundaries between clusters, and therefore measures linear separability of the point cloud, in an unsupervised manner. Prior knowledge was fed to the algorithm by setting the number of components to 𝑁 and by initializing the centroids on the positions of the primordia in a model without any noise:

∀𝑝 ∈ ⟦0, 𝑁⟦, { 𝜃 𝑝 0 ̂= 𝑝 ⋅ 𝛼 𝑟 𝑝 0 ̂= 𝑅 ⋅ 𝑒 𝛽𝑝
After convergence, the algorithm returns 𝑁 centroid points {(𝑟 𝑝 ̂, 𝜃 𝑝 ̂) ∶ 𝑝 ∈ ⟦0, 𝑁⟦} that we use to construct the Voronoi diagram. This actually defines a predictor for primordium rank by looking inside which cell of the diagram lies a given point (Supplementary Figure 1b):

𝑝̂𝑖 = argmin 𝑝∈⟦0,𝑁⟦ ‖(𝑟 𝑖 , 𝜃 𝑖 ) -(𝑟 𝑝 ̂, 𝜃 𝑝 ̂)‖
Finally, we estimate the Voronoi separability 𝑣 of our cloud of primordia points by computing the accuracy of the primordium rank prediction (noting 𝛿 𝑖 𝑗 the Kronecker delta on integers):

𝑣 = 1 𝑁 * 𝐾 ∑ 𝛿 𝑝 ̂𝑖 𝑝 𝑖 𝑁 * 𝐾 𝑖=1
If 𝑣 equals to 1, it means that the primordium points group into perfectly identifiable clusters. This can be interpreted as the fact that 𝐾 randomly sampled individuals can be superimposed perfectly, once they have been centered and aligned on their primordium which is the closest to the ℙ 0 stage. This will obviously be the case (as long as the motion coefficient 𝛽 remains realistic, typically <1) if no noise is introduced into the system. If 𝜎 𝛼 = 𝜎 𝑇 = 0, then all individuals proceed from the same regular exact pattern, and the only variability will be the one of the instant of sampling 𝜏 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 𝜎 𝛼 between 0° and 20° and 𝜎 𝑇 between 0 and 2 plastochrones, in a first time with the values 𝑅 = 30𝜇𝑚, 𝛼 = 137.51°, 𝛽 = 0.23 corresponding to actual measured data (Extended Data Figure 2d-e). As expected, increasing the angular variability creates more elongated clusters (Supplementary Figure 1c) that still appear separated. Yet, the Voronoi separation introduces confusion between neighboring primordia, more specifically making ℙ 𝑝 and ℙ 𝑝+3 overlap (Supplementary Figure 1f). Interestingly, when we increase the plastochrone variability (Supplementary Figure 1d), the confusion concerns rather ℙ 𝑝 and ℙ 𝑝+5 (Supplementary 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% (Supplementary Figure 1e).

The landscape of separability in the 𝜎 𝛼 × 𝜎 𝑇 parameter space gives an insight on the effects of variability on a population of individuals (Extended Data Figure 2f-g). With no surprise, primordia points appear to be less and less identifiable as azimuthal or plastochrone 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 (Extended Data Figure 2f, red contour). Interestingly, the azimuthal variability we measure on our observed SAMs is close to the upper bound of 𝜎 𝛼 on this optimal subspace.

If we fix the value of 𝜎 𝛼 to 6.7° as measured in our data (Extended Data Figure 2f, black vertical line), then the maximal plastochrone variability that is allowed for the separability to remain at 100% is close to 0.3. 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 plastochrone variability greater than 0.3, which would translate into an uncertainty on organ initiation times of nearly 3 hours. We have therefore demonstrated that the SAM achieves, at least, this level of rhythmic precision (though it might even be more precise) at the scale of a population (𝑁=21). This consequently validates our first order assumption that all considered SAMs are in a steady regime of development with the same plastochrone duration, and that the whole set of individual primordia of a given rank forms a homogeneous population in terms of developmental age.

Another interesting feature evidenced by this analysis is the influence of the motion speed of primordia separability. We explored the 𝛽 × 𝜎 𝑇 parameter space by fixing the value of 𝜎 𝛼 to the observed one, and varying the motion coefficient 𝛽 between 0 and 0.4 (Supplementary Figure 1h). It appears that lowering the speed reduces the maximal possible value of 𝜎 𝛼 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 plastochrone variability. For instance a value of 𝜎 𝑇 = 1 that translates into a separability of 95% when 𝛽 = 0.1 suddenly drops to a separability of only 50% if 𝛽 is increased up to 0.4. However increasing motion speed does not affect this much the maximal value of 𝜎 𝑇 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.

changes the tolerance of the system to azimuthal and plastochrone variability. High rhythmic precision is always required to achieve seamless superposition. Red contour indicates 100% separability, white contours every lower 5%. Black line indicates experimental value of speed coefficient.

Nuclei marker image quantification and aligned L1 signal maps

Confocal microscopy images saved as CZI files are processed using a complex computational pipeline that involves several image analysis, computational geometry and data manipulation steps, performed sequentially as in Supplementary Figure 2a.

Image reading: CZI files produced through the ZEN software 4 of the LSM-710 microscope are opened using a Python script 5 and split into independent channels that are saved separately as INR image files [START_REF] Douady | Phyllotaxis as a Dynamical Self Organizing Process Part II: The Spontaneous Formation of a Periodicity and the Coexistence of Spiral and Whorled Patterns[END_REF] . This operation preserves all the information contained in the raw image format. In the specific case of acquisitions for which both pPIN1:PIN1-GFP and pRPS5a:DII-VENUS are imaged, the close emission wavelengths causes the PIN1-stained cell membranes to appear in the DII nuclei images. In that only case, the PIN1 signal intensity is subtracted from the DII image channel before saving the file.

ROI cropping:

We use the polygonal selection tool of the ImageJ software to manually define a region of interest in all the slices of every image, and doing so, digitally dissect the outermost organs to get meristem images with at most 6 visible organs. This binary mask is then applied on all the channels, and masked channels are saved in separate INR files.

Automatic nuclei detection:

For each image, we perform an automatic detection of nuclei points based on the masked pRPS5a:TagBFP channel only. It consists of a 3D image stack ℐ Tag defined on a regularly spaced voxel grid {ℐ Tag (𝑥, 𝑦, 𝑧) ∶ (𝑥, 𝑦, 𝑧) ∈ 𝐼} that has a potentially different voxel size 𝑣 on each dimension so that:

(𝑣 𝑥 , 𝑣 𝑦 , 𝑣 𝑧 ) 𝐼 = •
The image is converted into a 4D Gaussian scale space 7 where ℬ 𝜎 (𝑥, 𝑦, 𝑧)is the discrete ball of radius 𝜎 in the image voxel grid:

ℬ 𝜎 (𝑥, 𝑦, 𝑧) = {(𝑥′, 𝑦′, 𝑧′) ∈ 𝐼 ∶ ‖(𝑥 ′ , 𝑦 ′ , 𝑧 ′ ) -(𝑥, 𝑦, 𝑧)‖ ≤ 𝜎}
The detection results in a 3D point cloud 𝒩 where each detected nuclei 𝑛 ∈ 𝒩 is represented by a single position 𝑃(𝑛) = (𝑥, 𝑦, 𝑧) in physical coordinates. 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. The parameter testing led to the determination of the optimal values 𝐾 𝜎 = 3, 𝜎 𝑚𝑖𝑛 = 0.8𝜇𝑚, 𝜎 𝑚𝑎𝑥 = 1.4𝜇𝑚, ℐ 𝑚𝑖𝑛 = 3000 corresponding to an evaluated performance of 95.6% recall and 98.5% precision.

Nuclei signal quantification: Every detected nucleus is assigned one value of signal intensity per image channel. This value is obtained by computing a weighted average of the masked channel intensity ℐ 𝑆 around the position of the nucleus. The signal images showing some local subcellular noise, the raw voxel value might not be fully representative of the whole nucleus. We chose to use a distance-based Gaussian weight, of constant radius 𝜎 𝑆 = 2µ𝑚 for all channels, to account for as much as possible of the signal information inside the nuclei, for which the typical measured diameter is ∼ 5𝜇𝑚. Practically, we first filter the signal image by a Gaussian kernel of radius 𝜎 𝑆 and read the values at the voxel positions of all detected nuclei: ∀𝑛 ∈ 𝒩 , 𝑆(𝑛) = ℐ 𝑆 𝜎 𝑆 (𝑃(𝑛))

Cell layer estimation: For the purpose of the analysis, we want to discriminate between the first layer of cells (L1) 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 as a 3D triangular mesh based on the channel. The image 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 8 on a resampled version of the image. This mesh undergoes a phase of triangle decimation 9 and isotropic remeshing 10 to obtain a surface composed of roughly 50000 regular faces. After normal estimation, only the largest connected component of triangles facing towards the upper side of the meristem is kept and used as the estimated meristem surface.

For each vertex 𝑣 of this surface mesh ℳ, we look for the closest detected nucleus by computing the distance to all of the points {𝑃(𝑛) ∶ 𝑛 ∈ 𝒩}. The nuclei that we label as "L1" are the set of nuclei points that are closer any other to at least one vertex of the surface:

∀𝑛 ∈ 𝒩 , 𝑛 * ∈ 𝐿 1 ⇔ ∃𝑣 ∈ ℳ ∶ 𝑛 * = argmin 𝑛∈𝒩 (‖𝑃(𝑛) -𝑃(𝑣)‖)
Rigid time registration: The previous steps were performed individually on each frame of the timelapse acquisitions. 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 12 applied once again on the pRPS5a:TagBFP channel of the consecutive images. This produces 𝒯 -1 isometry matrices in homogeneous coordinates 𝑅 𝑡 𝑖 ←𝑡 𝑖+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 these matrices to transform the images into the reference frame of the first frame of the sequence (𝑡 0 = 0ℎ). The transformed images are saved as separate INR files.

Registered nuclei points:

We also apply the obtained transforms to the positions of the nuclei detected at time 𝑡 𝑖 to obtain their registered positions 𝑃 0 in the reference frame of the first frame of the sequence:

∀𝑛 ∈ 𝒩(𝑡 𝑖 ), 𝑃 0 (𝑛) = ( ∏ 𝑅 𝑡 𝑗 ←𝑡 𝑗+1 0 𝑗=𝑖-1

) 𝑃(𝑛)

Non-linear time registration: In a second time, we want to estimate the local deformation of the 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 rigidly registered pRPS5a:TagBFP images, using again the block matching framework. This approach has proven to be efficient on plant tissues in the case of limited deformations 13,14 which is definitely the case for the meristematic tissue we are considering with a 4ℎ to 5ℎ interval between frames.

We do not use directly the resulting registered images but rather the vector field estimated in one of the two possible ways 𝑈 ⃗ ⃗ 𝑡 𝑖 →𝑡 𝑖+1 0 (respectively 𝑈 ⃗ ⃗ 𝑡 𝑖 ←𝑡 𝑖+1 0

) that stores 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 (respectively from the next frame to the current one). These two vector fields are saved as threedimensional INR files.

Nuclei cellular motion:

We estimate cellular motion between two consecutive sequence frames taken at times 𝑡 𝑖 and 𝑡 𝑖+1 in the forward direction using the first layer of registered nuclei points 𝑃 𝑡 𝑖 We deduce speed vectors measuring the local forward speed of cellular motion by dividing the motion vector by the time interval:

𝑣 𝑡 𝑖 →𝑡 𝑖+1 0 (𝑛) = 1 𝑡 𝑖+1 -𝑡 𝑖 𝑢 ⃗ 𝑡 𝑖 →𝑡 𝑖+1 0 (𝑛)
Conversely, similar cellular vectors are computed in the backward direction the 𝑈 ⃗ ⃗ 𝑡 𝑖 ←𝑡 𝑖+1 0 vector fields.

2D maps of epidermal signal:

We use a generic tool to take a signal defined on a discrete set of points and make it continuous in space. It basically computes a weighted average of signal values at every possible location in space. The weighting function we use is a parametric sigmoid density function 𝜂 𝑅,𝑘 of the distance 𝑟 to a point, which relies on two parameters, an extent parameter 𝑅 and a sharpness 𝑘 so that 𝜂 𝑅,𝑘 (0) ∼ 1, 𝜂 𝑅,𝑘 (𝑅) = 1 2 ⁄ and 𝜂̇𝑅 ,𝑘 (𝑅) =-𝑘 (Supplementary Figure 3a) and takes the form:

𝜂 𝑅,𝑘 (𝑟) = 1 2 - 1 2 tanh(𝑘 • (𝑟 -𝑅))
The continuous signal is then defined at every point 𝑝 of space based on a point cloud {𝑃(𝑛) ∶ 𝑛 ∈ 𝒩} and the associated signal values {𝑆(𝑛) ∶ 𝑛 ∈ 𝒩} as:

∀𝑝 ∈ ℝ dim(𝑃) , 𝑆 𝑅,𝑘 ⏟ 𝒩 (𝑝) = ( ∑ 𝜂 𝑅,𝑘 (‖𝑝 -𝑃(𝑛)‖) 𝑛∈𝒩 ) -1 ∑ 𝜂 𝑅,𝑘 (‖𝑝 -𝑃(𝑛)‖) 𝑛∈𝒩 • 𝑆(𝑛)
Note that, for any point 𝑝 where ∑ 𝜂 𝑅,𝑘 (‖𝑝 -𝑃(𝑛)‖) 𝑛∈𝑁 = 0, the signal map is not defined. To make the map outlines closer to their actual support, we will even consider that 𝑆 𝑅,𝑘 ̂(𝑝) is defined only if

∑ 𝜂 𝑅,𝑘 (‖𝑝 -𝑃(𝑛)‖) 𝑛∈𝑁 ≥ 1 2
⁄ . This constraint is equivalent to consider the implicit surface (or curve) obtained with the point cloud as generator and 𝜂 𝑅,𝑘 as potential as the boundary of the definition domain of 𝑆 𝑅,𝑘 ̂.

In the rest of the analysis, we will mostly compute 2D maps of epidermal signal, using the XY projection 𝑃 Π of a point cloud 𝑃 (𝑃 Π = {(𝑥, 𝑦) ∶ (𝑥, 𝑦, 𝑧) ∈ 𝑃}) and considering only the nuclei labeled as 𝐿1 (Supplementary Figure 3c-e). For instance the projected epidermal map at the time 𝑡 𝑖 of a sequence registered in the reference frame of the first image of the sequence will be computed as:

∀𝑝 ∈ ℝ 2 , 𝑆 𝑅,𝑘 0,Π ⏟ 𝐿 1 (𝑡 𝑖 ) (𝑝) = ( ∑ 𝜂 𝑅,𝑘 (‖𝑝 -𝑃 𝑡 𝑖 0,Π (𝑛)‖) 𝑛∈𝐿 1 (𝑡 𝑖 ) ) -1
∑ 𝜂 𝑅,𝑘 (‖𝑝 -𝑃 𝑡 𝑖 0,Π (𝑛)‖) 𝑛∈𝐿 1 (𝑡 𝑖 )

• 𝑆(𝑛)

To determine the best parameter values for the function 𝜂 𝑅,𝑘 we ran an extensive parameter exploration and measured the error made by mapping the signal to retain the values that yield the minimal error. This is measured by the average relative error between the actual signal value of a nucleus and the value, at the projected nuclei position 𝑃 0,Π (𝑛), of the 2D map computed with all nuclei but the considered one:

ℇ 𝑆 (𝐿 1 ) = 1 |𝐿 1 | ∑ | 𝑆 𝑅,𝑘 0,Π ⏟ 𝐿 1 ∖𝑛 (𝑃 0,Π (𝑛)) -𝑆(𝑛)| 𝑆(𝑛) 𝑛∈𝐿 1
The optimal values are the ones that minimize this error over the whole available set of SAM sequences 𝒮 = {{𝓈 𝑗 (𝑡 𝑖 ) ∶ 𝑖 = 0 ⋯ 𝒯 𝑗 } ∶ 𝑗 = 1 ⋯ |𝒮|} for the qDII signal, as it is our main focus in this work. The values we obtain, shown in Supplementary Figure 3b, are 𝑅 * = 7.5𝜇𝑚 and 𝑘 * = 0.55𝜇𝑚 -1 :

𝑅 * , 𝑘 * = argmin 𝑅,𝑘 (∑ ∑ ℇ qDII (𝐿 1 (𝓈 𝑗 , 𝑡 𝑖 )) 𝒯 𝑗 𝑖=0 |𝒮| 𝑗=1
)

From now on, we will consider that, if not explicitly mentioned otherwise, the maps referred to are 2D projected maps computed on all the epidermal nuclei using the optimal parameters. To simplify the notations we will then write:

𝑆 𝓈 𝑗 ,𝑡 𝑖 ̂= 𝑆 𝑅 * ,𝑘 * Π ⏟ 𝐿 1 (𝓈 𝑗 ,𝑡 𝑖 )
Sequence SAM reference frame determination: A common 3D cylindrical reference frame for the SAM sequences can be described by landmarking a set of key geometrical features for each meristem:  the unitary vector 𝕒 ⃗ ⃗ of the main vertical axis of the shoot apex  the position 𝕔 of the apex center in the central zone (CZ) of the meristematic dome  the unitary radial vector 𝕡 ⃗ ⃗ 0 ∶ 𝕒 ⃗ ⃗ ⊥ 𝕡 ⃗ ⃗ 0 of the direction of ℙ 0 relatively to 𝕔  where ℙ 0 is the position of the center of the last initiated organ primordium  the orientation 𝕠 ∈ {-1; 1} of the phyllotactic spiral (clockwise or counter-clockwise)

We estimate the center position 𝕔 on each sequence using the signal information from the pCLV3::mCHERRY channel quantified as {CLV3(𝑛) ∶ 𝑛 ∈ ⋃ 𝐿 1 (𝑡 𝑖 ) 𝒯 𝑖=0 }. In the map CLV3 0 ̂{𝑡 𝑖 ∶𝑖=0⋯𝒯}

the central zone appears as a wide isotropic peak of signal intensity (Supplementary Figure 3f). To extract the center and the radius of this peak, we threshold the map with a series of 𝐾 𝛾 values depending of the average signal intensity of the sequence CLV3 The use of a normalized CLV3 threshold made the values of 𝑟 𝕔 very homogeneous among the considered individuals around an average value of 28µm (Supplementary Figure 3g). Finally, the 𝑧 position of the CZ center is estimated using the local 𝑧 information:

𝑧 𝕔 = 𝑧 0 ̂{𝑡 𝑖 ∶𝑖=0⋯𝒯} (𝕔 Π )
This center can be placed before estimating the main vertical axis 𝕒 ⃗ ⃗ only because the tissues are imaged from above and we can make the approximation that 〈𝕒 ⃗ ⃗ , [0,0,1]〉 ∽ 1. In any case we try to position of this axis by applying a small rotation 𝑅 𝑥 (𝜓 𝑥 ) • 𝑅 𝑦 (𝜓 𝑦 ) around the x and y axes to the sequence registered positions centered on 𝕔 𝑃 𝕔 = 𝑃 0 -𝕔. This transformation leaves 𝕔 unchanged. We note:

𝑃 𝜓 𝑥 ,𝜓 𝑦 𝕔 = 𝑅 𝑥 (𝜓 𝑥 ) • 𝑅 𝑦 (𝜓 𝑦 ) • 𝑃 𝕔
Each rotation is associated with a main axis vector 𝕒 ⃗ ⃗ 𝜓 𝑥 ,𝜓 𝑦 = 𝑅 𝑥 (𝜓 𝑥 ) • 𝑅 𝑦 (𝜓 𝑦 ) • [0,0,1]. Given the radial symmetry of the meristematic dome, the optimal axis should be the one for which the function 𝑧(‖(𝑥, 𝑦)‖) 𝜓 𝑥 ,𝜓 𝑦

𝕔

shows less variability, at least up to the radial distance where the dome ends. We take this extent into account by applying a Gaussian weight of radius 𝜎 𝑟 = 20𝜇𝑚, function of the distance 𝑟 = ‖(𝑥, 𝑦)‖ of the rotated points to the center. We look for the optimal values of the rotation angles with the simple constraint that 〈𝕒 ⃗ ⃗ 𝜓 𝑥 ,𝜓 𝑦 , [0,0,1]〉 ≥ 0.95 by evaluating the error to the local mean 𝑧 of the rotated points:

𝜓 𝑥 𝕒 ⃗ ⃗ ,

𝕔

. The next step is to locate the 𝕡 ⃗ ⃗ 0 direction. This can be seen as a 2D problem if we consider the points 𝑃 𝕒 ⃗ ⃗ 𝕔 , for which the 𝑥𝑦 plane is orthogonal to 𝕒 ⃗ ⃗ .

If we consider a polar coordinate system centered on for this plane(𝑟, 𝜃), we are ultimately looking for angular coordinate 𝜃 ℙ 0 of the direction of the ℙ 0 primordium. In Supplementary Note 1, we detailed why ℙ 0 should correspond to a maximal auxin concentration in the meristem. Therefore, we look for the minimal value of the qDII ̂{𝑡 𝑖 ∶𝑖=0⋯𝒯}

𝕒 ⃗ ⃗ 𝕔 map, restraining the search to 𝑟 ∈ [𝜌 𝑚𝑖𝑛 𝑟 𝕔 ; 𝜌 𝑚𝑎𝑥 𝑟 𝕔 ] to avoid artifactual detections (we used 𝜌 𝑚𝑖𝑛 = 0.9 and 𝜌 𝑚𝑎𝑥 = 1.4 in the analysis)

𝜃 ℙ 0 = argmin 𝜃∈]-𝜋;𝜋] ( min 𝑟∈[𝜌 𝑚𝑖𝑛 𝑟 𝕔 ;𝜌 𝑚𝑎𝑥 𝑟 𝕔 ] qDII ̂{𝑡 𝑖 ∶𝑖=0⋯𝒯} 𝕒 ⃗ ⃗ 𝕔 (𝑟, 𝜃))
Finally, instead of estimating the orientation from the data, we chose to rely on manual expertise to determine visually the clockwise or counterclockwise orientations of the meristems. This produces a simple identity or reflection matrix, depending on the case:

𝑂 𝑧 (𝕠) = [ 1 0 0 0 𝕠 0 0 0 1 ]
Aligned L1 nuclei points: In the end, the determination of the SAM landmarks 𝕒 ⃗ ⃗ , 𝕔, 𝕡 ⃗ ⃗ 0 and 𝕠 allows to transform the sequence registered points 𝑃 0 into a common 3D reference frame in which we will be able to compare different individuals locally. We note 𝑃 * the positions of L1 nuclei points in this common reference frame:

∀𝑛 ∈ 𝐿 1 , 𝑃 * (𝑛) = 𝑅 𝑧 (𝜃 ℙ 0 ) • 𝑂 𝑧 (𝕠) • 𝑅 𝑥 (𝜓 𝑥 𝕒 ⃗ ⃗ ) • 𝑅 𝑦 (𝜓 𝑦 𝕒 ⃗ ⃗ ) • (𝑃 0 (𝑛) -𝕔)
Aligned L1 image slices: We also use the obtained transform to register the original images into the common reference frame by applying the same transform expressed in homogeneous coordinates:

[ 𝑅 𝑧 (𝜃 ℙ 0 ) • 𝑂 𝑧 (𝕠) • 𝑅 𝑥 (𝜓 𝑥 𝕒 ⃗ ⃗ ) • 𝑅 𝑦 (𝜓 𝑦 𝕒 ⃗ ⃗ ) 𝑇 𝕔 0 1 ]
This registration is performed using the blockmatching computational library 12 and creates a registered image ℐ 𝑆 * expressed over a voxel grid 𝐼 * that is centered on 0 in x and y and slightly shifted in z so that:

𝐼 * = 𝐼 -( 𝐾 𝑥 2 , 𝐾 𝑦 2 , 𝐾 𝑧 8 ) • (𝑣 𝑥 , 𝑣 𝑦 , 𝑣 𝑧 )
We used the 2D aligned 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 Figure 2a-c, Figure 3a-b, Figure 4a and Figure 5a-h. More specifically, we use the 𝑧 * coordinate of aligned 𝐿 1 nuclei points to compute the 2D map 𝑧̂ * over the centered grid coordinates 𝐼 * , and produce a 2D image defined the same xy grid as the original one but keeping the voxel value of one z slice per pixel:

ℐ 𝑆 * Π ⏟ 𝐿 1 = {ℐ 𝑆 * (𝑥, 𝑦, ⌊ 𝑧̂ * (𝑥, 𝑦) 𝑣 𝑧 ⌋ • 𝑣 𝑧 ) : (𝑥, 𝑦, . ) ∈ 𝐼 * }
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.

Organ primordia 2D detection:

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 try to use some a priori knowledge on regular phyllotaxis to detect the positions of the ranked organ primordia in the meristem, namely ℙ 0 , ℙ 1 , ℙ 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 qDII ̂ * map, but also that soon after organ initiation, an auxin depletion area is formed, creating a local maximum in qDII ̂ * . In that respect, our primordia detection procedure consists first in detecting extremal points in this map and in labeling them in a second time by organ primordium rank.

To locate extremal points in 2D, we are not only interested in absolute extremality (namely peaks and troughs of the map) but also for points that are extremal in a given direction, and we will therefore detect ridges and valleys to construct a landscape transform of the qDII ̂ * map. Peaks, which are maximal in any direction, will appear as convergence points of ridges, throughs as convergence points of valleys, and we will also retrieve saddle points as convergence points of a ridge and a valley.

Classically, extremal points are detected as zero-crossings of the gradient of a function. In our case, we discriminate between positive zero-crossings ∅ + and negative zero-crossings ∅ -that would correspond respectively to local maxima and local minima when applied on a gradient. For 1D scalar function 𝑓 ∶ ℝ → ℝ, zero-crossings are defined as follows:

∅ + 𝑓(𝑥) = (𝑓(𝑥) > 0) ∧ (𝑓(𝑥 + 𝑑𝑥) < 0) ∅ -𝑓(𝑥) = (𝑓(𝑥) < 0) ∧ (𝑓(𝑥 + 𝑑𝑥) > 0)
Now, if the scalar function is a 2D function 𝑓 ∶ ℝ 2 → ℝ, we can define a 2D zero-crossing, checking if the function is locally changing sign along each dimension, for instance:

∅ + ⃗⃗⃗⃗⃗ 𝑓(𝑥, 𝑦) = [ ∅ + 𝑓(⋅, 𝑦)(𝑥) ∅ + 𝑓(𝑥,⋅)(𝑦) ]
The map gradient is indeed a 2D function, but is not scalar since the gradient itself will be a 2D vector.

To cope with this, we introduce the notion of (here positive) zero-crossing of a function 𝑓 ⃗⃗⃗ ∶ ℝ 2 → ℝ 2 , along the direction given by the angle 𝜔, written ∅ 𝜔 + 𝑓 , measuring if the scalar function restricting 𝑓 ⃗⃗⃗ to the given direction is changing sign when we follow this direction:

∅ 𝜔 + 𝑓 (𝑥, 𝑦) = 〈∅ + ⃗⃗⃗⃗⃗ 〈𝑓 , 𝑒 𝜔 ⃗⃗⃗⃗ 〉(𝑥, 𝑦) , 𝑒 𝜔 ⃗⃗⃗⃗ 〉 where 𝑒 𝜔 ⃗⃗⃗⃗ = [ cos 𝜔 sin 𝜔 ]
Using this local directional zero-crossing information, we can compute ∅ + 𝑓 as a scalar measure of the proportion of directions in which the vector field is consistently changing sign, by integrating the zerocrossing of 𝑓 ⃗⃗⃗ over all possible directions:

∅ + 𝑓 (𝑥, 𝑦) = 1 2𝜋 ∫ ∅ 𝜔 + 𝑓 (𝑥, 𝑦) 𝜋 -𝜋

𝑑𝜔

By applying this integrated zero-crossing operator, we can compute the "landscape transform" of a map (Supplementary Figure 3h), consisting of a local ridge estimation Λ + and a local valley estimation Λ -, and obtain a way to quantify in detail local extremality on the qDII ̂ * map:

Λ + qDII ̂ * (𝑥, 𝑦) = ∅ + ∇ ⃗ ⃗ qDII ̂ * (𝑥, 𝑦) Λ -qDII ̂ * (𝑥, 𝑦) = ∅ -∇ ⃗ ⃗ qDII ̂ * (𝑥, 𝑦)
We first look for saddle points in order to eliminate them and disconnect the otherwise continuous networks of ridges and valleys. Saddle points  belong to both a valley and a ridge and are detected using their product, then valleys  and ridges  can be found outside the saddle locations, using different thresholds in both cases:

Λ  qDII ̂ * (𝑥, 𝑦) = √ Λ + qDII ̂ * (𝑥, 𝑦) • Λ -qDII ̂ * (𝑥, 𝑦) ≥ 𝑧 𝑚𝑖𝑛 Λ  qDII ̂ * (𝑥, 𝑦) = (Λ + qDII ̂ * (𝑥, 𝑦) ≥ 𝜆 𝑚𝑖𝑛 ) ∧ (¬Λ  qDII ̂ * (𝑥, 𝑦)) Λ  qDII ̂ * (𝑥, 𝑦) = (Λ -qDII ̂ * (𝑥, 𝑦) ≥ 𝜆 𝑚𝑖𝑛 ) ∧ (¬Λ  qDII ̂ * (𝑥, 𝑦))
We run this binary labeling with the values 𝑧 𝑚𝑖𝑛 = 0.01 and 𝜆 𝑚𝑖𝑛 = 0.03 and consider the connected components of each binary map as the set of extremal features of each kind, resulting in a set of regions (Supplementary Figure 3h), denoted for instance Λ 𝑙  for the connected valley labeled 𝑙. To sum up the information, each element 𝑙 of those three label sets written ,  and  is associated with a single spatial position (𝑟 𝑙 , 𝜃 𝑙 ), a measure of extremality 𝜆 𝑙 , a value of signal qDII ̂ * 𝑙 and an area 𝐴 𝑙 as such:

∀𝑙 ∈ ,
𝐴 𝑙 = ∬ 𝑑𝑟𝑑𝜃 Λ 𝑙  ∀𝑙 ∈ , { (𝑟 𝑙 , 𝜃 𝑙 ) = argmax (𝑟,𝜃)∈Λ 𝑙  (Λ + qDII ̂ * (𝑟, 𝜃) • Λ -qDII ̂ * (𝑟, 𝜃)) 𝜆 𝑙 = max (𝑟,𝜃)∈Λ 𝑙  (4 √ Λ + qDII ̂ * (𝑥, 𝑦) • Λ -qDII ̂ * (𝑥, 𝑦)) qDII ̂ * 𝑙 = qDII ̂ * (𝑟 𝑙 , 𝜃 𝑙 ) 𝐴 𝑙 = ∬ 𝑑𝑟𝑑𝜃 Λ 𝑙 
In a second time, we need to identify which of these extremal points are the ones corresponding to the youngest organ primordia. For this we use several assumptions from prior biological knowledge on the organization of the shoot apical meristem and the auxin distribution at the level of the primordium, namely:

(i) Organ primordia are marked by a local maximum of auxin concentration (ii) Primordia are organized in a regular spiral of divergence angle 𝛼 * = 2𝜋/𝜑 2 ~137.5° (iii) In this spiral, organ primordia can be ranked by their distance to the center (iv) No organ can initiate inside the central zone (CZ) of the meristem (v) The auxin depletion zone of a primordium is located between its maximum and the CZ If we consider the organ primordium of rank 𝑖 hereafter noted ℙ 𝑖 , we look for two extremal points, a maximum of auxin (minimal point of qDII ̂ * ) 𝑃(ℙ 𝑖  ) = (𝑟 ℙ 𝑖  , 𝜃 ℙ 𝑖  ) that actually defines the position of the primordium (i) and a minimum of auxin (maximal or saddle point of qDII ̂ * ) 𝑃(

ℙ 𝑖  ) = (𝑟 ℙ 𝑖  , 𝜃 ℙ 𝑖  )
which should follow the following rules:

(i) ∃ℙ 𝑖 ⟹ ∃𝑙 ∈  ∶ (𝑟 ℙ 𝑖  , 𝜃 ℙ 𝑖  ) = (𝑟 𝑙 , 𝜃 𝑙 ) (ii) 𝜃 ℙ 𝑖  ~𝜃ℙ 𝑖 * = 𝑖 • 𝛼 * (iii) 𝑟 ℙ 𝑖-1  ≤ 𝑟 ℙ 𝑖  ≤ 𝑟 ℙ 𝑖+1  (iv) 𝑟 ℙ 𝑖  ≥ 𝑟 ℂℤ (v) 𝑟 ℙ 𝑖  < 𝑟 ℙ 𝑖  , 𝜃 ℙ 𝑖  ~𝜃ℙ 𝑖 
To take these assumptions into account, we perform our identification of ℙ 𝑖  and ℙ 𝑖  among the detected extremal points  and  by associating to each of them a positive score 𝛾 ℙ 𝑖 (𝑙) relatively to ℙ 𝑖 , measuring how good a candidate it would be for ℙ 𝑖  (for points in ) or for ℙ 𝑖  (for points in  ∪ ). The idea behind the formula is to provide a result that complies as much as possible with the rules mentioned here above. The score is the product of an intrinsic, primordium independent, score 𝛾(𝑙) and of a spatial term 𝛾 ℙ 𝑖 (𝑙) = 𝛾 𝜃 ℙ 𝑖 (𝑙)

𝑖 -𝑑  ) • 𝑟 𝕔 ) 2 , 1) •(𝑟 𝑙 < 𝑟 ℙ 𝑖  ) if 𝑙 ∈  ∪ 
For each considered primordium ℙ 𝑖 , we first identify the point corresponding to a maximal value of auxin, provided there exists one with a score value above a 𝛾 min threshold:

ℙ 𝑖  = argmax {𝑙∈∶ 𝛾(𝑙)• 𝛾 ℙ 𝑖 (𝑙)>𝛾 min } (𝛾(𝑙) • 𝛾 ℙ 𝑖 (𝑙))
Then this detected point is used to compute the scores of minimal points of auxin by the means of 𝑟 ℙ 𝑖  (set to +∞ of ℙ 𝑖  is not defined) and we identify the minimal point associated with the primordium, provided it exists:

ℙ 𝑖  = argmax {𝑙∈∪∶ 𝛾(𝑙)• 𝛾 ℙ 𝑖 (𝑙)>𝛾 min } (𝛾(𝑙) • 𝛾 ℙ 𝑖 (𝑙))
Each identified primordium extremal point is then associated with a spatial position (𝑟 ℙ 𝑖  , 𝜃 ℙ 𝑖  ) or 3i) which can be used to compute estimates of all the spatial signals (starting with qDII ̂ * ) which allows us to track signals at the level of organ primordia.

(𝑟 ℙ 𝑖  , 𝜃 ℙ 𝑖  ) (Supplementary Figure
On all the meristems presented in the article, we preformed a manual correction of the primordium assignment of auxin extremal points, to make sure that we recover information from biologically meaningful locations. This manual labeling allowed us to perform an evaluation of the detection method, leading to an average 75.4% of correct detection (using a Jaccard index) rising to 87.3% when considering only primordia stages between ℙ -2 and ℙ 3 (Supplementary Figure 3j). Finally, using the fact that all 2D signal maps from different acquisitions are registered into the same reference spatial frame, we reconstruct the average map over a population of meristems by generalizing this formula over the whole available set of time-lapse sequences 𝒮 = {𝓈 𝑗 ∶ 𝑗 = 1 ⋯ |𝒮|} to combine all the quantitative spatio-temporal information into one dynamic map reflecting the canonical behavior of the system: This computation results in dynamic 2D maps as those presented for Auxin signal in Supplementary Video 1. These were obtained with parameter values 𝑅 𝜏 = 0.1, 𝑘 𝜏 = 2, 𝑝 min = -3, 𝑝 max = 5.

Developmental Time Estimation

𝑆 * ̂(𝑟,
are the ones used throughout the analyses. c-e. Using the L1 nuclei detected in the confocal image (c.) and their quantified signal values projected in 2D (d.) we compute a 2D map of epidermal signal (e.) in this case qDII. f. Using such a map computed on the quantified CLV3 levels, we are able to locate precisely the center of the meristem and estimate the extent of the CZ. g. The extents of CZ estimated using the CLV3 maps show a very limited variability around the 28µm value among the observed individuals (N=21 SAMs). h. Based on qDII maps a detection of extremal features (ridges, valleys and saddles) is performed and extremal points are extracted. i. Assuming a minimal regularity in position w/r to CZ and divergence angles close to 137.5°, the extremal points are given a score relatively to each considered primordium rank, and a maximum of one auxin maximal point and one auxin minimal point is assigned to each primordium. j. Comparison of automatically detected auxin extremal points with expertized ones demonstrate a very accurate detection between ranks P0 and P3, with a decreased performance when the features are less well defined (no absolute minimum before P0, several maxima after P4). Color indicates the rank of primordia. Filled color indicates accurate detection, light color indicates correct detection but inaccurate location, dark grey indicates false negatives, light grey indicates false positives.

Cell layer estimation from segmented images:

To be able to automatically determine to which layer a given cell belongs to, we use the topology of the tissue, notably the "background region" 𝒮 1 . Indeed, the biological definition of the epidermis, also called layer-1 or L1, is to be in contact with the outside world. We therefore label as 𝐿 1 any cell 𝑐 such that its region 𝒮 𝑐 can be considered a neighbor to 𝒮 1 .

Neighborhood relationship is defined through the notion of surface of contact between two regions in he segmented image. In a first time, we consider the 6-connectivity to define neighborhoods at voxelscale. However, to be robust to potential segmentation errors, we estimate the area of all surfaces of contact representing walls between two cells 𝑐 and 𝑐 ′ , and consider as neighbors only those for which the area is greater than a given threshold 𝐴 min . We compute the area of contact 𝐴 𝑐,𝑐 ′ 𝒮 as the sum of surfel areas between pairs of 6-adjacent voxels so that one is labeled 𝑐 and the other 𝑐 ′ . An analysis of the distribution of wall areas (data not shown) led us to consider that a value of 𝐴 min = 5𝜇𝑚 2 would be suitable. We note 𝑁(𝑐) the neighbors of the cell 𝑐, namely the labels that verify the surface of contact condition:

𝑁(𝑐) = {𝑐 ′ ∈ 𝒮, 𝐴 𝑐,𝑐 ′ 𝒮 > 𝐴 min }
Therefore, we obtain 𝐿 1 = 𝑁 (1). Subsequently, we can define the cells belonging to the second layer (𝐿 2 ) as those in contact with the 𝐿 1 , also filtered by the same minimal wall area:

𝐿 2 = ⋃ 𝑁(𝑐) ∩ 𝑐∈𝐿 1 (𝐿 1 ∪ {1}) ̅̅̅̅̅̅̅̅̅̅̅̅̅ .

Reconstruction of L1 anticlinal walls:

In order to quantify the PIN signal intensity at the level of each individual cell membrane, we need to have a precise identification of the cell walls and a faithful 3D representation to describe their position and orientation. If the boundary between two cells can be extracted as a set of voxels in the segmented image 𝒮, it will generally be too sensitive to noise and to resolution to be used as such. We chose therefore to use a triangular mesh representation with a high resolution to represent accurately the cell walls.

To obtain such meshes, we apply the Marching Cubes algorithm 8 to each cell 𝑐 represented by its connected region 𝒮 𝑐 of identically labeled voxels in the segmented membrane image 𝒮. This produces a triangular mesh ℳ 𝑐 , generally closed (except on image borders) and with voxel-like resolution. This mesh represents the shape of the cell by a set of vertices with 3D coordinates 𝑃 𝑐 and a set of triangular faces linking those vertices 𝑇 𝑐 : ℳ 𝑐 = 〈𝑇 𝑐 , 𝑃 𝑐 〉.

Using Marching Cubes ensures us that, in an 8-voxel cube with only two labels 𝑐 and 𝑐′ (which typically occurs at the interface between two cells) the algorithm will create the same vertices whichever label is considered as 1 or 0. In other words, we know that two cells that are neighbors in the image (with a large enough surface of contact) will have common vertices in their mesh reconstructions ℳ 𝑐 and ℳ 𝑐′ . We use this property to construct the cell wall mesh 𝒲 𝑐,𝑐′ = 𝒲 𝑐′,𝑐 of the interface between 𝑐 and 𝑐 ′ as the restriction of one of the two meshes to the set of common vertex points 𝑃 𝑐 ⋂𝑃 𝑐′ : 𝒲 𝑐,𝑐′ = 〈{𝑡 ∈ 𝑇 min(𝑐,𝑐 ′ ) ∶ ∀𝑣 ∈ 𝑡, 𝑃(𝑣) ∈ 𝑃 𝑐 ⋂𝑃 𝑐′ } , 𝑃 𝑐 ⋂𝑃 𝑐′ 〉 Each wall mesh undergoes then a phase of triangle decimation 9 and isotropic remeshing 10 to obtain a regular surface so that the typical length of a triangle edge is close to 0.5µm, which is about the voxel characteristic dimension. On the triangular mesh, we estimate the normal vectors 𝑛 𝑐→𝑐 ′ 𝑣 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ at each vertex, ensuring they all point from 𝑐 to 𝑐′, and the area of each triangle that allows us to estimate the total area 𝐴 𝑐,𝑐 ′ of the interface between 𝑐 and 𝑐′. Note that the Marching Cube intersection, along with the decimation and smoothing, will produce a mesh that is smaller that the actual wall (the intersecting part does not extend to cell edges) and the wall area will be underestimated. On the other hand, voxelbased area estimation is known to be largely overestimating, which ultimately provides a way to have both lower and upper bound estimates of the wall areas.

Quantification of PIN signal at wall-level:

We consider that the PIN polarity 𝕡 𝑐→𝑐 ′ = -𝕡 𝑐 ′ →𝑐 of a given cell interface, i.e. whether the PIN efflux carriers orient the flow of auxin from 𝑐 to 𝑐′ or from 𝑐′ to 𝑐, is given by the differential of PIN concentration that exists between the plasma membranes of 𝑐 and 𝑐′ at their interface. Our way to access this information is through the difference of PIN signal intensity in the image on either side of the cell wall marked by the PI signal intensity around 𝒲 𝑐,𝑐′ . We designed a method that relies on the sampling of both PIN and PI image signals across the wall to make a decision on whether there exists a significant difference in PIN intensity levels, and therefore a polarity of the wall.

The vertices and normals of the triangular mesh 𝒲 𝑐,𝑐′ are used to generate a set of 3D cylinders in which we will sample the image signals (Supplementary Figure 4a). We chose to exclude the vertices located on the contour of the mesh as the normal estimation may be less robust when the vertex is not fully surrounded by faces, and only consider the inside of the mesh 𝒲 𝑐,𝑐′ ∘ . At each inner vertex 𝑣, the main axis of the corresponding cylinder 𝒞 𝑐,𝑐′ 𝑣 is given by the local normal vector to the mesh 𝑛 𝑐→𝑐 ′ 𝑣 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ and is characterized by its radius 𝑟 𝒞 and its extent on each side of the wall 𝑑 𝒞 (Supplementary Figure 4bd). In the following we will consider only the image voxels lying inside this cylinder, which, if we note(𝑥, 𝑦, 𝑧) Π the projection of a voxel (𝑥, 𝑦, 𝑧) on the main axis (𝑥, 𝑦, 𝑧) Π = 𝑃(𝑣) + 〈(𝑥, 𝑦, 𝑧) -𝑃(𝑣), 𝑛 𝑐→𝑐 ′ 𝑣 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ 〉 𝑛 𝑐→𝑐 ′ 𝑣 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , can be defined as:

𝐼 ∩ 𝒞 𝑐,𝑐 ′ 𝑣 = {(𝑥, 𝑦, 𝑧) ∈ 𝐼 ∶ ‖(𝑥, 𝑦, 𝑧) Π -𝑃(𝑣)‖ ≤ 𝑑 𝒞 ∧ ‖(𝑥, 𝑦, 𝑧) -(𝑥, 𝑦, 𝑧) Π ‖ ≤ 𝑟 𝒞 }
We project the signal values of these voxels on an 1-dimensional axis assigning them a signed abscissa 𝑑 𝑣 (𝑥, 𝑦, 𝑧) = 〈(𝑥, 𝑦, 𝑧) -𝑃(𝑣), 𝑛 𝑐→𝑐 ′ 𝑣 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ 〉, which allows us to consider the spatial distribution of signal intensities along the cylinder axis (Supplementary Figure 4c). On this axis, the position 𝑑 𝑣 = 0 corresponds to the intersection of the cylinder 𝒞 𝑐,𝑐′ 𝑣 with the wall mesh 𝒲 𝑐,𝑐′ . This position might be shifted from the actual cell wall due to the stacking of residual errors from the consecutive processing steps (sequence segmentation artifacts, meshing simplifications, smoothing approximations). The idea is therefore to come back to the actual membrane-marker image signal ℐ 𝑃𝐼 (𝑥, 𝑦, 𝑧) to locate the position of the cell wall in the cylinder with a higher precision.

We achieve this precise location by estimating the abscissa of the mode of membrane signal intensity from the distribution {(𝑑 𝑣 , ℐ 𝑃𝐼 )(𝑥, 𝑦, 𝑧) ∶ (𝑥, 𝑦, 𝑧) ∈ 𝐼 ∩ 𝒞 𝑐,𝑐 ′ 𝑣 }. This is done by fitting a Gaussianshaped function, optimizing its parameters to minimize the squared error to the distribution (Supplementary Figure 4e). The function uses 4 parameters and its form is given by:

PI ̂(𝑑 𝑣 ) = PI 0 + PI 1 𝑒 - (𝑑 𝑣 -𝑑 0 PI ) 2 𝑑 1 PI 2
The least-squares estimation of the function notably gives us a value 𝑑 0 PI that marks the estimated abscissa of the PI maximum within the considered cylinder, which we will use as a reference defining the actual position of the cell wall. Then, the signals are quantified on either side of this abscissa, up to a fixed distance 𝑑 max , by computing the average voxel intensity (Supplementary Figure 4f). For instance, we determine the two values of PIN signal intensity as follows: 4g). In such case the polarity is given by the sign of the difference between the medians PIN 𝑐→𝑐 ′ and PIN 𝑐 ′ →𝑐 of the two distributions:

PIN 𝑐→𝑐 ′ 𝑣 = 1 |𝐼 ∩ 𝒞 ̇𝑐→𝑐 ′ 𝑣 | ∑ ℐ PIN (𝑥,
𝕡 𝑐→𝑐 ′ = { 1 if (p-value < 0.05) ∧ (PIN 𝑐→𝑐 ′ > PIN 𝑐 ′ →𝑐 ) -1 if (p-value < 0.05) ∧ (PIN 𝑐→𝑐 ′ < PIN 𝑐 ′ →𝑐 ) 0 otherwise
Note that, once again, the definition is symmetrical ensuring that 𝕡 𝑐→𝑐 ′ = -𝕡 𝑐 ′ →𝑐 . In any case, we also define the wall-level intensity of PIN signal as the average of the signal medians of either side:

PIN 𝑐,𝑐 ′ = PIN 𝑐 ′ ,𝑐 = 1 2 (PIN 𝑐→𝑐 ′ + PIN 𝑐 ′ →𝑐 )
Finally, each cell wall is also characterized by its barycenter 𝑃 𝑐,𝑐 ′ and by a single normal vector computed as the average of the inner wall normals:

𝑛 𝑐→𝑐 ′ ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ = -𝑛 𝑐 ′ →𝑐 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ = 1 |𝒲 𝑐,𝑐′ ∘ | ∑ 𝑛 𝑐→𝑐 ′ 𝑣 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ 𝑣∈𝒲 𝑐,𝑐′ ∘

Computation of cell polarity vectors:

The local wall polarity needs to be integrated to the cell level if we want to describe the local directionality of PIN carriers inside the tissue. To do this, we define the polarity of a cell 𝑐 as a 3D vector PIN 𝑐 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ that takes into account all the anticlinal cell walls surrounding the considered cell. Each wall contributes to the resulting vector proportionally to its area, and if polarized, adds up a directional flux parallel to its normal vector with an intensity equal to the difference of PIN intensities on its sides:

PIN 𝑐 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ = 1 ∑ 𝐴 𝑐,𝑐 ′ 𝑐 ′ ∈𝒩(𝑐)∩𝐿 1 ∑ 𝐴 𝑐,𝑐 ′ 𝑐 ′ ∈𝒩(𝑐)∩𝐿 1 (|PIN 𝑐→𝑐 ′ -PIN 𝑐 ′ →𝑐 |𝕡 𝑐→𝑐 ′ 𝑛 𝑐→𝑐 ′ ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ )
Alignment and map generation: To make it possible for the data computed on membrane images to be mapped onto the data from the nuclei channels of the same acquisition, we need to transform the quantified information into the same reference frame.

References Supplementary Method 3: Extrapolated cell motion in the developmental continuum

Throughout this work, we strongly relied on the spatio-temporal periodicity of phyllotactic systems, that we demonstrated to be a valid assumption for our considered SAMs. Indeed, the high angular precision and limited plastochrone variability make the observed systems close to a steady developing regular phyllotactic system with a divergence angle 𝛼 = 137.5 ± 6.7° and a plastochrone 𝑇 = 12 ± 2h (Extended Data Figure 2d-g, Supplementary Method 1).

Notably, we considered that, in the 2D cylindrical reference frame centered on the CZ of the shoot apical meristem, the dynamics of any quantifiable signal 𝑆 must follow the properties of a spatiotemporally periodic function of spatial period [0, -𝛼] and temporal period 𝑇: ∀𝑟 ∈ ℝ + , ∀𝜃 ∈ [-𝜋, 𝜋], ∀𝑡 ∈ 𝒯, 𝑆(𝑟, 𝜃, (𝑡 + 𝑇)) = 𝑆(𝑟, (𝜃 + 𝛼)[2𝜋], 𝑡) This helped us consider signal dynamics on durations that largely overpass the observation range of 10 to 14 hours, by applying successive rotations of the meristems to simulate the passing of time. Notably, an aligned SAM observed at 𝑡 rotated of 137.5° clockwise has been shown to be the best next frame to the same aligned SAM observed at 𝑡 =10h (Figure 1g, Extended Data Figure 2h). More generally, any primordium of stage 𝑝 visible at time 𝑡 can be use to infer information at time 𝑡 + 𝑝𝑇. This means that we were able to reconstruct trajectories of signals at a given position (𝑟, 𝜃) over a duration of up to 9 plastochrones (~100 hours) from observations spanning only one, but with 9 visible primordia (P-3 to P5). This was achieved only by interpolating rotated aligned information (Supplementary Movie 1).

Unfortunately, this global reconstruction heuristic could work only while we were looking at the same location in space, where the spatio-temporal property holds. To some extent, it can be generalized to robust primordia landmarks, such as auxin maxima, that we assume to be unique while moving in the course of primordium development. If they can be identified, and associated with a primordium of stage 𝑝, then they can be positioned on the same developmental axis at a time 𝑡 + 𝑝𝑇 to reconstruct a developmental history at the level of this time-tracked landmark.

However, the moment we are interested in cellular processes, such as the dynamics of transcriptional response to auxin for instance, the reconstructed long-term trajectories cannot be used to draw relevant conclusions, as they reflect the dynamics either at fixed coordinates or at non-cell-specific landmark points. It is therefore necessary to find a way to access temporal cell-level information. Individual cells can be tracked in time-lapse sequences, either manually or automatically, which could be use to obtain signal trajectories over 10 to 14 hours (Figure 2g). But to achieve the mentioned 100-hour reconstruction, spatio-temporal periodicity has to be used at some point.

Extrapolated tissue area tracking

We assume the cells in the central and peripheral zones (CZ and PZ respectively) of the SAM to have essentially an outward radial motion that accelerates as cells exit the central zone. This has been confirmed by the cellular motion vectors estimated from vector fields of image deformation (Supplementary Method 2) in which the azimuthal component is in average close to 0, with a limited amplitude compared to the radial component (Supplementary Figure 5a). We note 𝑣 * ⃗⃗⃗⃗ (𝑛) the speed of the 𝑛 in the 2D polar SAM reference frame:

∀𝑛 ∈ 𝐿 1 , 𝑣 * ⃗⃗⃗⃗ (𝑛) = 𝑣 𝑟 (𝑛)𝑟 * ⃗⃗⃗ (𝑛) + 𝑣 𝜃 (𝑛)𝜃 * ⃗⃗⃗⃗ (𝑛)

where 𝑟 * ⃗⃗⃗ (𝑛) and 𝜃 * ⃗⃗⃗⃗ (𝑛) are the local normal and tangential unitary vectors at 𝑃 * (𝑛). In the following, we will then assume a pure radial motion of cells in the L1 of the SAM, i.e. that:

∀𝑛 ∈ 𝐿 1 , 𝑣 𝜃 (𝑛)~0

We compute the local radial cellular speed as a 2D continuous map on the L1 (Supplementary Figure 4b-c), using the same parameters as for the other cellular signals: We perform this progression until we compute the last position 𝑃(𝑡 𝑛 ) for which motion can not be estimated (as there is no next image it compute image deformation) However, we can still extrapolate the lastly computed motion to reach one plastochrone, and estimate 𝑃(𝑡 0 + 𝑇) with:

𝑟(𝑡 0 + 𝑇) = 𝑟(𝑡 𝑛-1 ) + (𝑇 -(𝑡 𝑛-1 -𝑡 0 )) ⋅ 𝑣 𝑟 * ̂(𝑟(𝑡 𝑛-1 ), 𝜃 0 , 𝑡 𝑛-1

To proceed further, we would like to progress in time and estimate the cell position at 𝑡 1 + 𝑇. This is where we use the spatio-temporal periodicity property to derive that: We are therefore able to estimate this next position from observations at 𝑡 0 simply by rotating the radial speed map (Supplementary Figure 5d). Then iteratively it is possible to go further in time to 𝑡 𝑛 + 𝑇, extrapolate motion to 𝑡 0 + 2𝑇, apply again spatio-temporal periodicity to reach 𝑡 1 + 2𝑇, and so on. These are valid for positive integer values of 𝑝 until reaching the maximal extent of 𝑃 in the considered data (i.e. when the map used to estimate local radial speed becomes ill-defined), which determines a value 𝑝 max . This defines a radial trajectory in the 2D space that reflects the local motion of cells along several plastochrones (Supplementary Figure 5e). Note that a rigorously identical approach can be used to go backwards in time with negative values of 𝑝 until 𝑡 0 + 𝑝 min 𝑇, using motion vectors computed using inverse image deformation (Supplementary Method 2).

In the end, from an initial position (𝑟 0 , 𝜃 0 , ) we obtain a discrete radial trajectory:

{(𝑟(𝑡 𝑖 + 𝑝𝑇), 𝜃 0 ) ∶ 𝑖 ∈ ⟦0, 𝑁⟧, 𝑝 ∈ ⟦𝑝 min , 𝑝 max ⟧}

To monitor a cellular process over long time courses, the objective would be to estimate the value of the signal 𝑆 along this spatio-temporal trajectory, i.e:

{𝑆 * ̂(𝑟(𝑡 𝑖 + 𝑝𝑇), 𝜃 0 , 𝑡 𝑖 + 𝑝𝑇) ∶ 𝑖 ∈ ⟦0, 𝑁⟧, 𝑝 ∈ ⟦𝑝 min , 𝑝 max ⟧} Using the spatio-temporal periodicity property of 𝑆, this translates into:

{𝑆 * ̂(𝑟(𝑡 𝑖 + 𝑝𝑇), 𝜃 0 + 𝑝𝛼, 𝑡 𝑖 ) ∶ 𝑖 ∈ ⟦0, 𝑁⟧, 𝑝 ∈ ⟦𝑝 min , 𝑝 max ⟧}

In other terms, we have defined a series of spatial locations in a time-series of acquisitions such that the sequence of signal values at these locations on the same meristem estimates the cell-level trajectory of the considered signal. If we represent them on a time-series of meristems, we define tissue areas that can be tracked, first in time then in space, to reconstruct the average behavior of a group of cells over time (Supplementary Figure 5f-h). This is the approach we use to reconstruct celllevel auxin trajectories (Figure 2h) and to study the relationship between auxin input and transcriptional response in a consistent group of cells (Figure 4d, Figure 4f-h).

Figure 1 :

 1 Figure 1: Spatial auxin distribution in the SAM follows a precise and reiterative pattern.

Figure 2 .

 2 Figure 2. Auxin dynamics are not cellularly fixed

Figure 3 .

 3 Figure 3. Auxin polar transport and biosynthesis act synchronously to regulate patterning.

Figure 4 .

 4 Figure 4. Auxin and its transcriptional output show a complex non-linear relationship.

Figure 5 .

 5 Figure 5. Signal temporal integration regulates auxin output through epigenetic control.

Extended Data Figure 1 .

 1 Expression pattern of qDII and CLV3. a-c. Representative images of RPS5a:qDII. TagBFP channel (grays). d. Quantification of TagBFP nuclei intensity at the L1 cell layer of the meristem (N=42991 nuclei). e-g. DII-VENUS expression pattern. h. Quantification of DII-VENUS-N7 nuclei intensity at the L1 cell layer. i-k. Expression pattern of pCLV3:mCherry. i. Quantification of CLV3 nuclei intensity at the L1 cell layer (N=6003 nuclei). Scale bars = 20 µM. m. Comparison of CLV3 expression as a function of radial distance from the center in time. Each point represents a nuclei and regression curves for each time point distribution. Boxplot centers show median. Extended Data Figure 2. Precision of maxima positioning enables time extrapolation a. Superposition of non-registered SAM images and registered by rotation. N=69 meristem images. b. Superposition of the same images in (a) after rigid registration. Scale bars = 20 µM. c. Distributions of relative map errors comparing an individual qDII map to the average of the rest of the population. From left to right, the boxplots represent the errors obtained using the aligned image positions as in (b); the aligned image positions but only in the P0 domain; the aligned image positions compared with 90°-rotated individual maps; the raw image positions centered on the CLV3 domain; and the raw image positions as in (a). The alignment significantly reduces the error, not only in the P0 domain that is used for alignment but all over the SAM, making around 3 times less error that the worst possible alignment. d. Auxin maxima radial position as a function of developmental time, assuming a plastochrone time of 12h. Radial distance of an auxin maximum Pn at T= 0h lies between the previous maximum Pn-1 at T=10h and T=14h. Boxplot centers show median. e. Relative divergence angle between 2 consecutive maxima. Boxplot centers show median. f. Theoretical study shows that variability vastly affects the separability of primordia clusters. In a phyllotactic model using the observed initial distance and speed (from d), divergence angle and angular variability (from e) parameter values the apparent separability of primordia across individuals can only be explained by very limited plastochrone variability (<0.3). Seamless superposition of individuals proves that the SAM achieves a very high rhythmic precision at population-scale. Red contour indicates 100% separability, white contours every lower 5%. Black line indicates experimental value of azimuthal variability. g. Primordia points obtained with a phyllotactic model for different values of angular and rhythmicity variability. Color marks the rank of the primordium. Symbols indicate the position in the parameter space (f) h. Computation of error between a map at T=10h and rotated maps (-180° to +180°) at T=10h (red), T=5h (green) or T= 0h (blue). Error values equal to 0 denote identical maps. The lowest error values with rotated maps (blue points) were globally found between the 10h map and the 0h map close to a 137° rotation, validating that the extrapolation heuristic of placing P1 at T=0h right after P0 at T=10h is indeed the most consistent way to extend a sequence in time Extended Data Figure 3. Auxin depletion areas are locally defined and transient in time. a. Auxin minima appear from stage P1 as a saddle point isolating radially the auxin maximum from the central zone, and progressively builds up until its reduction to a single cell file in the organ boundary at stage P5. b. Auxin values of maxima (filled boxes) and minima (outline boxes) in time for each primordia. c. Minima (outline boxes) closely follow their associated maximum (filled boxes) but stop at a fixed distance to the center after stage P3. d. Auxin depletion radial extend from the maxima to the center of the SAM. The depletion zone progressively widens in stages P1 to P4 before getting split in two at stage P4. e. Auxin depletion angular extent from the auxin maxima. Boxplot centers show median. Extended Data Figure 4. Croissant-shape PIN1 localization does not always reflect actual cell polarities. a. Confocal projection of PIN1-GFP. Magnified regions of the SAM in (b) and (i) are indicated by white squares. Scale bar = 20 µM. b. Confocal projection of PIN-signal. The arrow indicates the polarity for cell 2 assigned by visual impression based on intensity and croissant-shape of the signal. Scale bar = 2 µM. c-f. Z-slices of image in (b) with PIN1-GFP (green) and PI (magenta) signals. Notably the PIN1 signal belongs to the membrane of cell 1 rather than cell 2 (y shift between signals). g. Average signal distribution for PIN1 (green) and PI (grey) in the anticlinal wall between cell 1 and 2. Each line represents the average signal of a region of the anticlinal cell wall between cells. h.

  Quasi-static PIN1 network have a global convergence towards the center of the SAM. a. Representative confocal projection of the L1 cell layer of qDII (yellow), PIN1-GFP (green) and CLV3 (magenta). Average (N= 4 SAMS) auxin map and vector fields (b) showing convergence at the auxin maxima (black areas). c. local PIN1 convergence map. Blue represents convergence while brown represents divergent vectors. d. PIN1 expression levels map. e. Circumferential auxin distribution at a 35 µM radial distance (line) with corresponding PIN1 vector fields (green arrows). a-e. T= 0 h. f-j. T= 5 h. k-o. T= 10 h. p-t. T= 14 h. Extended Data Figure 6. Expression patterns of YUCCA genes in the SAM a-k. Expression patterns of YUCCA genes in the SAM. Transcriptional reporter lines for the YUCCA auxin biosynthetic genes (orange) are indicated. Propidium iodide staining is shown in gray. Scale bars = 20 µM. Extended Data Figure 7. Precision of DR5 radial positioning and stage-dependent Pearson correlations a. Radial distribution of CLV3 (magenta circles) and DR5 (black circles) at different locations of the SAM (P-1 to P1). N= 21 SAM. Each point represent a nuclei. Vertical line indicates the mean for DR5 radial distribution. b-c. Pearson correlation coefficients computed on tracked cell-level values of auxin input levels and output. The considered variables are the levels at time T=0h, T=5h and T=10h, the sum of auxin input from T=0h to T=5h and from T=0h to T=10h, and the variations of DR5 levels from T=0h to T=5h and from T=0h to T=10h. At a global level (b), no significant positive correlations are found between input and output variables, even when considering the time difference present in the data, arguing against a simple effect of delay, even between input levels and output variations, excluding a simple integrated response. Looking locally at nuclei in P-1 developmental stage makes stronger correlation values appear, and clear positive values (c). This overall suggests the existence of an exogenous stage-dependent variable controlling the relationship between auxin concentration and transcriptional response. Extended Data Figure 8. Time integrated auxin response is control by chromatin remodeling a. Confidence intervals (shade) and regression (line) for DR5 expression upon different auxin treatment lengths. Curves represent DR5 expression patterns at different regions of the PZ (ordinate) before and after treatment. Mock N= 1025 nuclei, 30 min N= 952 nuclei, 120 min N= 1344 nuclei, 300 min N= 2022 nuclei. b. Confidence Intervals and regression for DR5 expression of control (magenta) or TSA (green) treated meristems. DMSO N= 4184 nuclei, TSA N= 4141 nuclei.

Property 8 :

 8 Ordering on a regular phyllotaxis: In a system 𝒮 with a steady development of plastochrone 𝑇 and a regular phyllotaxis of divergence angle 𝛼, if 𝛼 is not a simple fraction of 2𝜋, i.e if: ∄𝑎, 𝑏 ∈ ⟦1, |𝒮|⟧, 𝛼 = 𝑎 𝑏 2𝜋 then the primordia angles 𝜃 𝑝 are sufficient to know the strict total order on 𝒮: ∀𝑝, 𝑞 ∈ 𝒮, (∃𝑡 ∈ 𝒯, ∃𝑘 ∈ ℤ, 𝜃 𝑞 (𝑡) -𝜃 𝑝 (𝑡) = 𝑘𝛼 * [2𝜋]) ⇒ 𝑞 = 𝑝 + 𝑘In this case, we say that 𝒮 has a clear regular phyllotaxis.

  𝑛) ∶ 𝑛 ∈ 𝐿 1 (𝑡 𝑖 )} and the transformation 𝑈 ⃗ ⃗ 𝑡 𝑖 →𝑡 𝑖+1 0 that maps the current frame into the next one by a vector field. Each nucleus can then be assigned a local motion vector 𝑢 ⃗ 𝑡 𝑖 →𝑡 𝑖+1 0 (𝑛) = 𝑈 ⃗ ⃗ 𝑡 𝑖 →𝑡 𝑖+1 0 (𝑃 𝑡 𝑖 0 (𝑛))

{(

  𝑟 𝑙 , 𝜃 𝑙 ) = argmin (𝑟,𝜃)∈Λ 𝑙  (qDII ̂ * (𝑟, 𝜃)) 𝜆 𝑙 = max (𝑟,𝜃)∈Λ 𝑙  (Λ -qDII ̂ * (𝑟, 𝜃)) + Λ -qDII ̂ * (𝑟 𝑙 , 𝜃 𝑙 ) qDII ̂ * 𝑙 = min (𝑟,𝜃)∈Λ 𝑙  (qDII ̂*(𝑟, 𝜃)) 𝜃)∈Λ 𝑙  (Λ + qDII ̂ * (𝑟, 𝜃)) + Λ + qDII ̂ * (𝑟 𝑙 , 𝜃 𝑙 ) qDII ̂ * 𝑙 = max (𝑟,𝜃)∈Λ 𝑙  (𝑞𝐷𝐼𝐼 ̂ * (𝑟, 𝜃))

1 ∑ 1 •

 11 We also use the same method to compute an aligned vectorial map of PIN polarities, this time based on the estimation of cell polarity vectors: PIN ⃗⃗⃗⃗⃗⃗⃗ * ̂(𝑝) = ( ∑ 𝜂 𝑅 * ,𝑘 * (‖𝑝 -𝑃 * (𝑐)‖) 𝜂 𝑅 * ,𝑘 * (‖𝑝 -𝑃 * (𝑐)‖) 𝑐∈𝐿 PIN 𝑐 * ⃗⃗⃗⃗⃗⃗⃗⃗⃗ In particular, this is the data we use to quantify the convergence of PIN directions at tissue scale by applying the divergence operator to the vector field resulting from the map computation PIN ⃗⃗⃗⃗⃗⃗⃗ * ̂(𝑥, 𝑦) = (PIN 𝑥 * ̂, PIN 𝑦 * ̂)(𝑥, 𝑦). This operation gives us a scalar map where negative values correspond to areas of local convergence of the PIN directions: div PIN ⃗⃗⃗⃗⃗⃗⃗ * ̂(𝑥, 𝑦) = ∇ • PIN ⃗⃗⃗⃗⃗⃗⃗ * scalar map is computed from the estimation of the vector field on a regular 2D grid of coordinates by applying a 1-D Gaussian derivative filter of standard deviation 𝜎 div = 𝑅 * /2 to each component of the vectorial map in the adequate dimension. We use this computed map to approximate the values of div PIN ⃗⃗⃗⃗⃗⃗⃗ * ̂ at punctual locations such as nuclei points 𝑃(𝑛) or primordia extremal points 𝑃( ℙ 𝑖 ), 𝑃( ℙ 𝑖 ) or 𝑃( ℙ 𝑖 ) by averaging the values at the 4 closest grid coordinates. Indeed, it would be very complex to estimate this quantity on a discrete set of points from the PIN 𝑐 * ⃗⃗⃗⃗⃗⃗⃗⃗⃗ vectors, on which the divergence operator would be non-trivial.

1 ∑ 1 ⋅

 11 ,𝑘 * (‖(𝑟, 𝜃) -𝑃 * (𝑛)‖) 𝑛∈𝐿 𝜂 𝑅 * ,𝑘 * (‖(𝑟, 𝜃) -𝑃 * (𝑛)‖) 𝑛∈𝐿 𝑣 𝑟 (𝑛) Defined this way, local radial cellular speed is a tissue-level information, that is not attached to a cell but to a spatial position (𝑟, 𝜃). Therefore it has a spatio-temporal periodicity property and we can write:𝑣 𝑟 * ̂(𝑟, 𝜃, 𝑡 + 𝑇) = 𝑣 𝑟 * ̂(𝑟, 𝜃 + 𝛼, 𝑡)We use this spatio-temporal periodicity property of local cellular speed to extrapolate cell motion over time from a series of acquisitions of SAMs at discrete times {𝑡 0 < ⋯ < 𝑡 𝑛 < 𝑇 }. Let us consider a cell with an initial position 𝑃(𝑡 0 ) = (𝑟(𝑡 0 ), 𝜃 0 ), setting 𝑟(𝑡 0 ) = 𝑟 0 . Using acquisitions at 𝑡 0 , it is possible to estimate 𝑣 𝑟 * ̂(𝑟 0 , 𝜃 0 , 𝑡 0 ), which we use to estimate 𝑃(𝑡 1 ) = (𝑟(𝑡 1 ), 𝜃 0 ) assuming a linear motion between 𝑡 0 and 𝑡 1 :𝑟(𝑡 1 ) = 𝑟 0 + (𝑡 1 -𝑡 0 ) ⋅ 𝑣 𝑟 * ̂(𝑟 0 , 𝜃 0 , 𝑡 0 )More generally, with observations at 𝑡 𝑖-1 , we derive 𝑃(𝑡 𝑖 ) with: 𝑟(𝑡 𝑖 ) = 𝑟(𝑡 𝑖-1 ) + (𝑡 𝑖 -𝑡 𝑖-1 ) ⋅ 𝑣 𝑟

𝑟(𝑡 1 +

 1 𝑇) = 𝑟(𝑡 0 + 𝑇) + (𝑡 1 -𝑡 0 ) ⋅ 𝑣 𝑟 * ̂(𝑟(𝑡 0 + 𝑇), 𝜃 0 , 𝑡 0 + 𝑇) = 𝑟(𝑡 0 + 𝑇) + (𝑡 1 -𝑡 0 ) ⋅ 𝑣 𝑟 * ̂(𝑟(𝑡 0 + 𝑇), 𝜃 0 + 𝛼, 𝑡 0 )

  Finally we obtain the two following general equations: ∀𝑖 ∈ ⟦1, 𝑛⟧, 𝑟(𝑡 𝑖 + 𝑝𝑇) = 𝑟(𝑡 𝑖-1 + 𝑝𝑇) + (𝑡 𝑖 -𝑡 𝑖-1 ) ⋅ 𝑣 𝑟 * ̂(𝑟(𝑡 𝑖-1 + 𝑝𝑇), 𝜃 0 + 𝑝𝛼, 𝑡 𝑖-1 ) 𝑟(𝑡 0 + (𝑝 + 1)𝑇) = 𝑟(𝑡 𝑛-1 + 𝑝𝑇) + (𝑇 -(𝑡 𝑛-1 -𝑡 0 )) ⋅ 𝑣 𝑟 * ̂(𝑟(𝑡 𝑛-1 + 𝑝𝑇), 𝜃 0 + 𝑝𝛼, 𝑡 𝑛-1 )

Extended Data Table 1. RNAseq expression of YUCCA genes at the SAM

	AGI	Name	TPM Col	TPM Col	TPM Col	TPM Col
			Rep1	Rep2	Rep3	mean
	AT5G11320 YUC4 7.56	10.31	12.17	10.02
	AT4G32540 YUC1 8.56	8.59	11.98	9.71
	AT5G25620 YUC6 1.50	2.37	1.94	1.94
	AT4G13260 YUC2 0.68	0.77	0.55	0.67
	AT1G04610 YUC3 0.05	0.08	0.15	0.09
	AT4G28720 YUC8 0.06	0.06	0.09	0.07

  At every time 𝑡 ∈ 𝒯, each primordium 𝑝 ∈ 𝒮 is characterized by its current state {𝜏 𝑝 (𝑡), 𝑥 𝑝 (𝑡), 𝑦 𝑝 (𝑡)} where: o 𝜏 𝑝 (𝑡) ∈ [𝜏 min , 𝜏 max ] ⊂ ℝ is called the developmental state of the primordium. o 𝑥 𝑝 (𝑡) = [𝑟 𝑝 (𝑡), 𝜃 𝑝 (𝑡), 𝑧 𝑝 (𝑡)] ∈ ℝ 3 is the spatial position of the primordium in a cylindrical coordinate system, the origin of which is called the center of the system. o 𝑦 𝑝 (𝑡) ∈ ℝ 𝑑 is a vector describing the physiological state of the primordium.  The developmental state 𝜏 𝑝 is a continuous strictly increasing function of time. Note that 𝜏 𝑝 is consequently a bijection between 𝒯 and 𝜏 𝑝 (𝒯).  When it exists, the time 𝑡 0

𝑝 ≝ 𝜏 𝑝 -1 (0) is called initiation time of the primordium 𝑝.

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

  In a system 𝒮 with a steady development of plastochrone 𝑇, we can extend the notation of the consecutiveness relationship in the strict total order of 𝒮 using the isomorphism 𝓈 𝑡 to identify the primordia by their relative developmental stages write: We say that a system 𝒮 with a steady development of plastochrone 𝑇 has a regular phyllotaxis of divergence angle 𝛼 if the constant parts of spatial positions of two consecutive primordia in the strict total order of 𝒮 only differ by a rotation of angle 𝛼 around the vertical axis:

	∀𝑝, 𝑞 ∈ 𝒮, ∃𝑘 ∈ ℤ, ∀𝑡 ∈ 𝒯 ∶ 𝓈 𝑡 (𝑞) = 𝓈 𝑡 (𝑝) + 𝑘	def ⇔ 𝑞 = 𝑝 + 𝑘
	0 𝛼 Definition 4: ∀𝑝, 𝑝 + 1 ∈ 𝒮, 𝑋 𝑝+1 = 𝑅 𝑧 (𝛼) ⋅ 𝑋 𝑝 = 𝑋 𝑝 + [
		0

∀𝑡 ∈ 𝒯, 𝜏 𝑝+1 (𝑡) = 𝜏 𝑝 (𝑡) + 1 Property 4: In a system 𝒮 with a steady development of plastochrone 𝑇, at any time 𝑡 ∈ 𝒯, the rounding function of the developmental state: 𝓈 𝑡 ∶ 𝒮 ⟶ ⟦𝓈 min (𝑡), 𝓈 max (𝑡)⟧ ⊂ ℤ 𝑝 ⟼ ⌊𝜏 𝑝 (𝑡) + 0.5⌋ is an isomorphism. We call 𝓈 𝑡 (𝑝) the developmental stage of primordium 𝑝 at time 𝑡. If 𝓈 𝑡 (𝑝) = 𝑘 ∈ ℤ we say that the primordium 𝑝 has the label ℙ 𝑘 at time 𝑡. Property 5: In a system 𝒮 with a steady development of plastochrone 𝑇, two consecutive primordia in the strict total order of 𝒮 always have consecutive developmental stages: ∀𝑝, 𝑝 + 1 ∈ 𝒮, ∀𝑡 ∈ 𝒯, 𝓈 𝑡 (𝑝 + 1) = 𝓈 𝑡 (𝑝) + 1 Definition 6: ] Property 7: Spatio-temporal periodicity: In a system 𝒮 with a steady development of plastochrone 𝑇 and a regular phyllotaxis of divergence angle 𝛼, the system verifies the following spatio-temporal periodicity: ∀𝑝, 𝑝 + 1 ∈ 𝒮, ∀𝑡 ∈ 𝒯, { 𝑥 𝑝+1 (𝑡) = 𝑅 𝑧 (𝛼) ⋅ 𝑥 𝑝 (𝑡 + 𝑇) 𝑦 𝑝+1 (𝑡) = 𝑦 𝑝 (𝑡 + 𝑇)

Problem 1: Assignation of developmental stages:

  Given a system 𝓢, observed at 𝒏 discrete time points {𝒕 𝟏 , … , 𝒕 𝒏 } in which, for every 𝒑 ∈ 𝓢 and for every 𝒕 𝒊 ∈ {𝒕 𝟏 , … , 𝒕 𝒏 }, only a partially observed state 𝒑 ̃(𝒕 𝒊 ) = {𝒙 ̃𝒑(𝒕 𝒊 ), 𝒚 ̃𝒑(𝒕 𝒊 )} is available; Find for every 𝑡 𝑖 ∈ {𝑡 1 , … , 𝑡 𝑛 } an estimated developmental stage function 𝓈̂𝒕 𝒊 that verifies:

	∀𝑝, 𝑝 + 1 ∈ 𝒮, 𝓈̂𝒕 𝒊 (𝑝 + 1) = 𝓈̂𝒕 𝒊 (𝑝) + 1

  by filtering it successively with 𝐾 𝜎 3D isotropic Gaussian filters 𝒢 of increasing radii, varying geometrically between 𝜎 𝑚𝑖𝑛 and 𝜎 𝑚𝑎𝑥 : Nuclei are detected as local 4D maxima in this scale space representation noted {ℐ Tag 𝜎 = ℐ Tag * 𝒢 𝜎 ∶ 𝜎 ∈ ⟦𝜎 𝑚𝑖𝑛 ; 𝜎 𝑚𝑎𝑥 ⟧ 𝑔 𝐾 𝜎 }, with the limit that only the maxima where the response Tag 𝜎 is higher than a threshold ℐ 𝑚𝑖𝑛 are retained. A 4D point (𝑥, 𝑦, 𝑧, 𝜎) is then considered a local maximum if and only if:

	⟦𝜎 𝑚𝑖𝑛 ; 𝜎 𝑚𝑎𝑥 ⟧ 𝐾 𝜎 = {𝜎 𝑚𝑖𝑛 ( 𝐾 ∀(𝑥 ′ , 𝑦 ′ , 𝑧 ′ ) ∈ ℬ 𝜎 (𝑥, 𝑦, 𝑧), ∀𝜎 ′ ∈ ⟦𝜎 𝑚𝑖𝑛 ; 𝜎 𝑚𝑎𝑥 ⟧ 𝑔 𝜎 𝑚𝑎𝑥 𝜎 𝑚𝑖𝑛 ) 𝑖 𝜎 𝑔 𝐾 𝜎 , ℐ Tag 𝜎 (𝑥, 𝑦, 𝑧) ≥ ℐ Tag 𝜎′ (𝑥 ′ , 𝑦 ′ , 𝑧 ′ )

𝜎 ∶ 𝑖 𝜎 = 0 ⋯ 𝐾 𝜎 } = 𝑒 ⟦ln(𝜎 𝑚𝑖𝑛 ), ln(𝜎 𝑚𝑎𝑥 )⟧ 𝐾 𝜎

  • ⟦𝛾 𝑚𝑖𝑛 ; 𝛾 𝑚𝑎𝑥 ⟧ 𝐾 𝛾 , and keep for each value the largest connected component (lcc) to compute its barycenter and area:

	ℂℤ 𝛾 = lcc(CLV3 0 ̂{𝑡 𝑖 ∶𝑖=0⋯𝒯} ≥ 𝛾CLV3)
	We retrieve the average barycenter and values of all such CZ domains obtained with 𝛾 𝑚𝑖𝑛 =
	1.2, 𝛾 𝑚𝑎𝑥 = 1.8 and 𝐾 𝛾 = 7:
	𝕔 Π = (𝑥; 𝑦) ∶ (𝑥; 𝑦) ∈ ℂℤ 𝛾 : 𝛾 ∈ ⟦𝛾 𝑚𝑖𝑛 ; 𝛾 𝑚𝑎𝑥 ⟧ 𝐾 𝛾
	𝑟 𝕔 = √ 𝐴(ℂℤ 𝛾 ) 𝜋	∶ 𝛾 ∈ ⟦𝛾 𝑚𝑖𝑛 ; 𝛾 𝑚𝑎𝑥 ⟧ 𝐾 𝛾

  • 𝛾 𝑟 ℙ 𝑖 (𝑙) that takes into account the organization of the primordia.

			(1 -qDII ̂ * 𝑙 ) • min (	1 2	+ 𝜆 𝑙 ,1) • (1 -	1 𝐴 𝑙	) • (𝑟 𝑙 > 𝑟 𝕔 2 ⁄ ) if 𝑙 ∈ 
		𝛾(𝑙) =	(qDII ̂ * 𝑙 ) • min ( 1 2	+ 𝜆 𝑙 ,1) • (1 -	1 𝐴 𝑙	) • (𝑟 𝑙 > 𝑟 𝕔 2 ⁄ )	if 𝑙 ∈ 
			{ (0.2 + qDII ̂ * 𝑙 ) • min ( 1 2	+ 𝜆 𝑙 ,1) • (𝑟 𝑙 > 𝑟 𝕔 2 ⁄ )	if 𝑙 ∈ 
			𝛾 𝜃 ℙ 𝑖 (𝑙) = cos (𝜋 • min (	|𝜃 𝑙 -𝜃 ℙ 𝑖 * | 𝜃 max ℙ 𝑖	, 1))
	𝛾 𝑟 ℙ 𝑖 (𝑙) =	min (( 𝑟 max ℙ 𝑖 • 𝑟 𝕔 𝑟 𝑙 min (( { (𝑟 max ℙ 𝑖 -𝑑  ) • 𝑟 𝕔 ) 2 , 1) • min (( 𝑟 min 𝑟 𝑙 ℙ 𝑖 • 𝑟 𝕔 𝑟 𝑙 ) 2 , 1) • min (( (𝑟 min 2 ) ℙ	, 1) 𝑟 𝑙	if 𝑙 ∈ 

:

  It is possible to position the different individual acquisitions from various time-lapse sequences onto the same developmental time axis. Making the assumption that all the considered meristems develop at a comparable rate of 2 new organs per day, corresponding to a plastochrone time of 12h, and that they can, at the first order, be considered synchronous (similarly labeled organ primordia being at comparable developmental stages) we propose a very simple developmental time indexation. For a sequence of 𝒯 acquisitions taken at {𝑡 0 = 0h, 𝑡 1 , ⋯ , 𝑡 𝒯-1 } we choose to assign to each acquisition at 𝑡 𝑖 a developmental time: Using the indexation of individual acquisitions in developmental time and the aligned L1 signal maps, we reconstruct a large-scale picture of the typical behavior of signals in the meristem by computing averages among individuals and interpolations over time.At the scale of one time-lapse sequence 𝓈 𝑗 , for which we have computed the 2D aligned signal maps {𝑆 𝓈 𝑗 ,𝑡 𝑖 * ̂∶ 𝑖 = 0 ⋯ 𝒯 𝑗 -1} we approximate the continuous time signal by diffusing the information in developmental time by the same method we used to diffuse signal in space when computing continuous maps. At any time position 𝜏, the estimated map is a weighted average of single time maps, where the weights are time-distance-based density coefficients computed using the 𝜂 𝑅,𝑘 function with a time radius 𝑅 𝜏 and a time slope 𝑘 𝜏 : However, this approximation only enables us to reconstruct the evolution of signals over a time range equivalent to the duration of acquisitions. To extrapolate further more and obtain a more complete view of primordium development, we use the spatio-temporal periodicity of the system to consider the information at different spatial locations as information at different temporal positions. More precisely, we use the fact that to look 𝑝 plastochrones further in time is equivalent to rotate the system of the angle 𝜃 ℙ 𝑝 corresponding to the direction of the primordium ℙ 𝑝 . Therefore we use the angles 𝜃 ℙ 𝑝  estimated with 2D primordium detection to apply rotations on the computed maps for primordia stages going from ℙ 𝑝 min to ℙ 𝑝 max and therefore compute a dynamic map that covers a much larger temporal range (the first formula becoming a particular case with 𝑝 min = 𝑝 max = 0):

	𝜏 𝑖 = 𝜂 𝑅 𝜏 ,𝑘 𝜏 (|𝜏 𝑖 -𝜏|) 12h 𝑡 𝑖 𝒯 𝑗 -1 1 ∑ 𝜂 𝑅 𝜏 ,𝑘 𝜏 (|𝜏 𝑖 -𝜏|) • 𝑆 𝓈 𝑗 ,𝑡 𝑖 * ̂ (𝑟, 𝜃) 𝑖=0 ∑ 𝒯 𝑗 -1 ∑ 𝑖=0 ∑ ̂(𝑟, 𝜃, 𝜏) = * L1 Dynamic Signal Maps: 𝑆 𝓈 𝑗 𝑆 𝓈 𝑗 * 𝜂 𝑅 𝜏 ,𝑘 𝜏 (|(𝜏 𝑖 + 𝑝) -𝜏|) • 𝑆 𝓈 𝑗 ,𝑡 𝑖 * ̂ (𝑟, 𝜃 + 𝜃 ℙ 𝑝  ) 𝒯 𝑗 -1 𝑝 max 𝑝=𝑝 min 𝑖=0 ̂(𝑟, 𝜃, 𝜏) = ∑ ∑ 𝜂 𝑅 𝜏 ,𝑘 𝜏 (|(𝜏 𝑖 + 𝑝) -𝜏|) 𝒯 𝑗 -1 𝑝 max 𝑝=𝑝 min 𝑖=0

  PI < 𝑑 max }. Note that the definition is symmetrical, meaning that the result would be exactly the same if 𝑐 and 𝑐′ were to be permutated.By performing this two-sided estimation on every cylinder defined by the wall triangular mesh, we end up with two parallel signal distributions {PIN 𝑐→𝑐 ′ We test statistically whether these distributions can be seen as significantly different by an ANOVA test, and decide that a polarity exists when the test gives a p-value < 0.05 (Supplementary Figure

						𝑣	∶ 𝑣 ∈ 𝒲 𝑐,𝑐 ′ ∘ } and {PIN 𝑐 ′ →𝑐 𝑣	∶ 𝑣 ∈ 𝒲 𝑐,𝑐 ′ ∘ }.
							𝑦, 𝑧)
						(𝑥,𝑦,𝑧)∈𝐼∩𝒞 ̇𝑐→𝑐 ′ 𝑣
		PIN 𝑐 ′ →𝑐 𝑣	=	1 |𝐼 ∩ 𝒞 ̇𝑐′ →𝑐 𝑣	|	∑ (𝑥,𝑦,𝑧)∈𝐼∩𝒞 ̇𝑐′ →𝑐 𝑣	ℐ PIN (𝑥, 𝑦, 𝑧)
	where we note 𝐼 ∩ 𝒞 ̇𝑐→𝑐 ′ 𝑣	= {(𝑥, 𝑦, 𝑧) ∈ 𝐼 ∩ 𝒞 𝑐,𝑐 ′ 𝑣 ∶ -𝑑 max < 𝑣 (𝑥, 𝑦, 𝑧) -𝑑 0 PI < 0 } and 𝐼 ∩ 𝒞 ̇𝑐′ →𝑐 𝑣	=
	{(𝑥, 𝑦, 𝑧) ∈ 𝐼 ∩ 𝒞 𝑐,𝑐 ′ 𝑣 ∶ 0 < 𝑑 𝑣 (𝑥, 𝑦, 𝑧) -𝑑 0

  The processes of Rigid time registration and Sequence SAM reference frame determination ultimately define a single rigid transformation matrix 𝑇 * per acquired time point, that transforms the physical image coordinates into coordinates in the common reference SAM frame. Therefore we simply need to apply this transform to the wall centers ∀𝑐 ∈ 𝐿 1 , ∀𝑐 ′ ∈ 𝒩(𝑐) ∩ 𝐿 1 , 𝑃 𝑐,𝑐 ′ 𝑃(𝑐), and to apply the rotation component 𝑅 * of this matrix to the cell polarity vectors ∀𝑐 ∈ 𝐿 1 ,

	PIN 𝑐 * ⃗⃗⃗⃗⃗⃗⃗⃗⃗ = 𝑅 * PIN 𝑐 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ .
	Using the epidermal map formalism, introduced in 2D maps of epidermal signal, and from these
	transformed coordinates, we can then compute an aligned continuous map of PIN intensity based of
	the wall-level quantification of wall signal and the estimation of wall areas:
			-1
	PIN * ̂(𝑝) =	∑	𝐴 𝑐,𝑐 ′ 𝜂 𝑅 * ,𝑘 * (‖𝑝 -𝑃 𝑐,𝑐 ′ * ‖)	∑	𝐴 𝑐,𝑐 ′ 𝜂 𝑅 * ,𝑘 * (‖𝑝 -𝑃 𝑐,𝑐 ′ * ‖)	• PIN 𝑐,𝑐 ′
		𝑐∈𝐿 1		𝑐∈𝐿 1
	(	𝑐 ′ ∈𝒩(𝑐)∩𝐿 1	)	𝑐 ′ ∈𝒩(𝑐)∩𝐿 1

* = 𝑇 * 𝑃 𝑐,𝑐 ′, to the cell barycenters ∀𝑐 ∈ 𝐿 1 , 𝑃 * (𝑐) = 𝑇 *

Acknowledgements

We thank Fabrice Besnard and the members of the SIGNAL team for insightful discussions; Antoine Larrieu for helping with RNA-Seq analysis; Hélène Robert-Boisivon for the YUC transcriptional lines; We acknowledge the contribution of SFR Biosciences (UMS3444/CNRS, US8/Inserm, ENS de Lyon, UCBL) facilities PLATIM, for assistance with microscopy. This work was supported by Human Frontier Science Program organization (HFSP) Grant RPG0054-2013 to T.V. and C.G., ANR-12-BSV6-0005 grant (AuxiFlo) to T.V. and DFG FOR2581 to J.L.

Data availability

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

Code availability

Generic quantitative image and geometry analysis algorithms are provided in Python libraries timagetk, cellcomplex and tissue_nukem_3d (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 (https://gitlab.inria.fr/gcerutti/sam_spaghetti.git).

All other custom source codes and analysis scripts are available from the corresponding authors upon request.

Author Contribution

Membrane marker image and PIN1 polarity quantification

Automatic segmentation of membrane images: In order to be able to quantify membrane-localized signal, we need a segmentation of the tissue at the cellular level. Such automatic cell segmentation procedures are often limited by their capacity to detect the right number of "seeds" prior to the segmentation of the image according to membrane localized signal.

In our case, the presence of a constitutive nuclei targeted signal (pRPS5a::TagBFP), allows to compute the nuclei coordinates in the image. It is thus possible to use these coordinates as seeds to initialize a watershed segmentation algorithm. The quality of the obtained segmentation, in terms of "correct number of cells detected", is then directly linked to the nuclei detection quality. Compared to parametric seed detection by methods such as local minima detection (by h-transform algorithm) followed by connected component labeling, the use of detected nuclei signal coordinates allows to reduce the over and under segmentation problems (data not shown). To summarize the pipeline used for this automatic membrane segmentation step, we performed:

1. An adaptative histogram equalization 15 for all z-slices of the membrane stacks to improve and normalize the contrast; 2. Isometric resampling to a voxelsize of (0.2, 0.2, 0.2µm), when original images are (0.2, 0.2, 0.5µm), to performs Gaussian smoothing and obtain smoother segmentation in Z; 3. Gaussian smoothing of the membrane intensity image, with 𝜎 = 0.2𝜇𝑚 to reduce noise in the image when performing watershed segmentation; 4. Create a "seed image" from the nuclei coordinates to initialize the watershed algorithm; 5. Run the seeded-watershed algorithm with isometric smoothed intensity image and seed image.

The steps 2 to 5 correspond to the pipeline described in 13 where more technical details can be found. No post-segmentation corrections where performed, no cell-fusion (in case of over-segmentation) or morphological corrections (median filters to smooth the walls). In the end we obtain a segmented image 𝒮 that assigns an integer label to every voxel of 𝐼. The cells of the tissue are represented by independent connected regions of voxels (so that the same label can not be assigned to voxels that are not part of the same connected component of 𝒮) . The background corresponds to a specific label, that is systematically set to 1 to ensure consistency between images. Each cell labeled 𝑐 > 1 is the then represented by a connected region 𝒮 𝑐 so that:

Cell barycenter extraction: To obtain the cells barycenter, we used the SciPy library 16 and the center of mass function that simply estimates the position 𝑃(𝑐) of the center of the cell labeled 𝑐 in the segmented image as the average of the voxel coordinates of the cell region 𝒮 𝑐 :