1. Developmental Biology
  2. Plant Biology
Download icon

Temporal integration of auxin information for the regulation of patterning

  1. Carlos S Galvan-Ampudia
  2. Guillaume Cerutti
  3. Jonathan Legrand
  4. Géraldine Brunoud
  5. Raquel Martin-Arevalillo
  6. Romain Azais
  7. Vincent Bayle
  8. Steven Moussu
  9. Christian Wenzl
  10. Yvon Jaillais
  11. Jan U Lohmann
  12. Christophe Godin  Is a corresponding author
  13. Teva Vernoux  Is a corresponding author
  1. Laboratoire Reproduction et Développement des Plantes, Univ Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Inria, France
  2. Department of Stem Cell Biology, Centre for Organismal Studies, Heidelberg University, Germany
Research Article
  • Cited 2
  • Views 2,810
  • Annotations
Cite this article as: eLife 2020;9:e55832 doi: 10.7554/eLife.55832

Abstract

Positional information is essential for coordinating the development of multicellular organisms. In plants, positional information provided by the hormone auxin regulates rhythmic organ production at the shoot apex, but the spatio-temporal dynamics of auxin gradients is unknown. We used quantitative imaging to demonstrate that auxin carries high-definition graded information not only in space but also in time. We show that, during organogenesis, temporal patterns of auxin arise from rhythmic centrifugal waves of high auxin travelling through the tissue faster than growth. We further demonstrate that temporal integration of auxin concentration is required to trigger the auxin-dependent transcription associated with organogenesis. This provides a mechanism to temporally differentiate sites of organ initiation and exemplifies how spatio-temporal positional information can be used to create rhythmicity.

eLife digest

Plants, like animals and many other multicellular organisms, control their body architecture by creating organized patterns of cells. These patterns are generally defined by signal molecules whose levels differ across the tissue and change over time. This tells the cells where they are located in the tissue and therefore helps them know what tasks to perform.

A plant hormone called auxin is one such signal molecule and it controls when and where plants produce new leaves and flowers. Over time, this process gives rise to the dashing arrangements of spiraling organs exhibited by many plant species. The leaves and flowers form from a relatively small group of cells at the tip of a growing stem known as the shoot apical meristem.

Auxin accumulates at precise locations within the shoot apical meristem before cells activate the genes required to make a new leaf or flower. However, the precise role of auxin in forming these new organs remained unclear because the tools to observe the process in enough detail were lacking.

Galvan-Ampudia, Cerutti et al. have now developed new microscopy and computational approaches to observe auxin in a small plant known as Arabidopsis thaliana. This showed that dozens of shoot apical meristems exhibited very similar patterns of auxin. Images taken over a period of several hours showed that the locations where auxin accumulated were not fixed on a group of cells but instead shifted away from the center of the shoot apical meristems faster than the tissue grew. This suggested the cells experience rapidly changing levels of auxin. Further experiments revealed that the cells needed to be exposed to a high level of auxin over time to activate genes required to form an organ. This mechanism sheds a new light on how auxin regulates when and where plants make new leaves and flowers. The tools developed by Galvan-Ampudia, Cerutti et al. could be used to study the role of auxin in other plant tissues, and to investigate how plants regulate the response to other plant hormones.

Introduction

Specification of differentiation patterns in multicellular organisms is regulated by gradients of biochemical signals providing positional information to cells (Rogers and Schier, 2011; Wolpert, 1969). In plants, graded distribution of the hormone auxin is not only essential for embryogenesis, but also for post-embryonic development, where it regulates the reiterative organogenesis characteristic of plants (Dubrovsky et al., 2008; Vanneste and Friml, 2009; Benková et al., 2003). 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; Figure 1A). In Arabidopsis thaliana, as in a majority of plants, organs are initiated sequentially in the SAM peripheral zone (PZ surrounding the CZ) at consecutive relative angles of close to 137°, either in a clockwise or anti-clockwise spiral (Figure 1AGalvan-Ampudia et al., 2016). SAM organ patterning or phyllotaxis has been extensively analyzed using mathematical models (Douady and Couder, 1996; Mitchison, 1977; Veen and Lindenmayer, 1977). A widely accepted model proposes that the time interval between organ initiations (the plastochron) and the spatial position of organ initiation emerge from the combined action of inhibitory fields emitted by pre-existing organs and the SAM center (Douady and Couder, 1996). Tissue growth then self-organizes organ patterning by moving organs away from the stem cells and leaving space for new ones.

Figure 1 with 3 supplements see all
Spatial auxin distribution in the SAM follows a precise reiterative pattern.

(A) SAM radial organization. The CZ (magenta) is surrounded by the PZ (cyan). Emerging flower primordia and flowers are colored in yellow. (B). Representative expression patterns of DII-VENUS-N7 (yellow) and pCLV3:mCherry transcriptional reporter line (magenta). Primordia are indicated by color and rank. (C). Auxin map (1-qDII, yellow to black) of (B). CLV3 expression (magenta) and radial extension (circle) are shown. Black arrows depict radial distance from the center and aligned angle. (D–F). Superposition of 21 aligned SAM images at time 0 hr (D), and 10 hr (E). (F) 137.5° clockwise rotation of (D) results in a quasi-identical image of (E). See Figure 1—figure supplement 2A for non-aligned image superposition. Scale bars = 20 µm. (G). Precision in auxin maxima positioning measured using angular position deviation (azimuthal deviation, left panel) and radial direction (right panel). Red lines indicate the average cellular distance. N = 21 meristems. Colors indicate primordium ranks (P-1 blue, P0 cyan, P1 green, P2 yellow, P3 orange) (H). Space can be used as a proxy for time, as a rotation of 1 divergence angle is equivalent to a translation of 1 plastochron in time.

Auxin is the main signal for positional information in phyllotactic patterning (Reinhardt et al., 2003a; Reinhardt et al., 2000). Auxin, has been proposed to be transported directionally toward incipient primordia where it activates a transcriptional response leading to organ specification (Benková et al., 2003; Reinhardt et al., 2003a; Heisler et al., 2005; Vernoux et al., 2000). PIN-FORMED1 (PIN1) belongs to a family of auxin efflux carriers whose polarity determines the direction of auxin fluxes (Benková et al., 2003; Gälweiler et al., 1998). PIN1 proteins are present throughout the SAM and regulate the spatio-temporal distribution of auxin cooperatively with other carriers (Reinhardt et al., 2003a; Bainbridge et al., 2008). Convergence of PIN1 carriers toward sites of organ initiation was proposed to control an accumulation of auxin that triggers organ initiation. This spatial organization of PIN1 polarities was also proposed to deplete auxin around organs, locally blocking initiation and thus establishing auxin-based inhibitory fields (Reinhardt et al., 2003a; Heisler et al., 2005; Vernoux et al., 2011; de Reuille et al., 2006; Stoma et al., 2008; Jonsson et al., 2006; Smith et al., 2006a). In addition, a reduced responsivity of the CZ to auxin has been demonstrated, providing an auxin-dependent mechanism for the inhibition of organogenesis in the CZ (Vernoux et al., 2011; de Reuille et al., 2006). Several models converge to suggest that together, these auxin-dependent regional cues determine new organ locations in the growing SAM.

The genetically-encoded biosensor DII-VENUS, a synthetic protein degraded directly upon sensing of auxin, recently allowed an unprecedented qualitative visualization of spatial auxin gradients in the SAM (Vernoux et al., 2011; Brunoud et al., 2012). However, quantification of the spatio-temporal dynamics of auxin is required to fully evaluate both experimental and theoretical understanding of the action of auxin in SAM patterning. This is all the more important given that the continuous helicoidal reorganization of auxin distribution in the growing SAM, suggests that auxin might convey complex positional information. Here, we used a quantitative imaging approach to question the nature of the auxin-dependent positional information. We further investigate how efflux and biosynthesis regulate the 4D dynamics of auxin, and explore how this information is processed in the SAM to generate rhythmic patterning.

Results

Spatio-temporal auxin distribution

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

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

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

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

Video 1
Auxin developmental continuum over nine plastochrons.

Auxin distribution dynamics in the SAM obtained from population averaging and temporal extrapolation. The developmental stage indicated at the top p=n corresponds to the area located on the right. Color code: yellow = low auxin, to black = high auxin.

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

Auxin information travels as centrifugal waves in the meristem.

(A–C) Representative projection of P1 nuclei showing DII-VENUS-N7 (yellow) and TagBFP nuclei (grey) intensity changes in time. Time tracked nuclei are marked by white and red circles showing rapid decrease or increase of auxin over 10 hr, respectively. The green circle is centered on the position of the auxin maxima at each time point. The magenta line indicates the limit of the CLV3 domain. Scale bars = 20 µm. (D) Average motion of maxima (colored lines) is faster than average cell motion (grey lines). The magenta line indicates the CLV3 domain border. N = 21 meristems. (E) Compared distributions of radial motion speeds of auxin maxima (color boxplots) estimated as the slope of a linear regression per individual. Individual nucleus radial speed (gray boxplots) at the location of the maxima. N = 21 meristems. (F–G) Individual cells experience different auxin histories. Tracked cells at different locations (F; colored circles) and corresponding auxin levels (ordinate) over time (0, 5, 10, 14 hr). Scale bar = 20 µm. (H) Cellular mean auxin trajectories as a function of radial distance. Each line represents an extrapolated cell-size sector moving accordingly to cellular radial motion by its Gaussian average trajectory in radial distance (abscissa) and auxin value (ordinate). The color indicates the developmental stages at a given radial distance (P-1 = blue, P0 = cyan, P1 = green, P2 = yellow).

Spatio-temporal control of auxin efflux and biosynthesis

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

Figure 3 with 3 supplements see all
Spatio-temporal organization of auxin fluxes and biosynthesis.

(A) Co-visualization of PIN1-GFP (green), DII-VENUS-N7 (yellow) and pCLV3:mCherry (magenta). Scale bar = 20 µm. Square shows P-1 sector. (B) Magnified P-1 region of (A) PIN1-GFP (green) and PI (magenta). (C) Computed PIN1 cell interface polarities of (B). Green arrows indicate polarities with a p-value<0.1, small arrows < 0.25 and dots > 0.25. (D–G). Image of PIN1-GFP (green) and cell wall (magenta) obtained using confocal (D) or super resolution (SRRF) microscopy (F) and respective PIN1 cell interface polarities (E,G). (H) Quantification of PIN1-GFP expression. N = 4 meristems. (I) PIN1 vector map (green arrows) organization correlated with auxin distribution (yellow to black). N = 4 meristems. (J) PIN1 polarity divergence index at auxin maxima (color filled boxplots) or auxin minima (white filled boxplots) positions during organ initiation. N = 4 meristems. (K) The YUC4 auxin biosynthesis limiting enzyme is specifically expressed in developing flowers. YUC4:GFP transcriptional reporter in yellow, cell wall (PI) staining in grey. Scale bars = 20 µm. (L) yuc1yuc4 mutant inflorescence and meristem morphological defects (inset). Scale bars are 10 mm and 20 µm (inset). (M) Schematic representation of the tissue-scale organization of auxin transport and biosynthesis in relation to auxin distribution.

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

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

The role of time in transcriptional responses to auxin

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

Figure 4 with 1 supplement see all
Auxin and its transcriptional output show a complex non-linear relationship.

(A–C). Time-lapse images of representative projections of DII-VENUS-N7 (yellow), TagBFP nuclei (grey) and pDR5:mTurquoise2 (cyan). Scale bars = 20 µm. (D). Quantified DR5 expression map (black to cyan). Colored sectors show the tissue areas where primordia are located (P-3 to P1). N = 21 meristems. (E). Principal Component Analysis (PCA) showing absence of correlation (orthogonality) between auxin and DR5 at the tissue scale. Colored ellipses show the consistent pattern associated with each primordium stage (from P-2 to P1 using the same colors as in (D)). (F). Auxin and DR5 non-linear relationship in primordia. Cells from P-3 to P1 are indicated with the color code used in (D). Lines represent the regression of auxin and DR5 medians in time. (G–H). Auxin (G) and DR5 (H) expression in primordia. Boxplots use the same color code for primordia as in (D). N = 21 meristems. 

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

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

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

Figure 5 with 2 supplements see all
Temporal integration of auxin concentration regulates transcription.

(A–I) Activation of the DR5 reporter with different concentrations of auxin and durations of treatments. pDR5:mTurquoise2 expression before auxin treatment (A–D) and 5 hr after the end of the auxin treatment: mock (E) or 1 mM IAA treatment for 30’ (F), 120’ (G) or 300’ (H). pDR5:mTurquoise (cyan), TagBFP driven by pRPS5a (gray) and pCLV3:mCherry (magenta) labelled nuclei are shown. Inset: DII-VENUS-N7 (yellow) from the same meristem. Quantification of DR5 expression in the PZ after auxin treatments. (I). Average DR5 response in the PZ with different auxin concentrations and treatment durations. Confidence intervals (shade) and regression (line) shows log(DR5) expression along the circumference of the PZ (aligned angle) of control (gray) or IAA (color) treated meristems. For simplicity only the angular position of primordia are indicated (in grey, presumptive positions). (J–M). Transcriptional response to auxin treatment of different durations in pid-14. pid-14 pDR5:3xVENUS SAM treated with IAA for 120’ (J,K) or 300’ (L,M) are shown. (N–Q). Transcriptional response to a 300’ auxin treatment in ett mutants. Control Col-0 pDR5:3xVENUS-N7 (N,O) and ett-22/arf3 pDR5:3xVENUS-N7 (P,Q) meristems treated with auxin for 300’.

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

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

Discussion

In a recent modeling study, a stochastic induction of organ initiation based on temporal integration of morphogenetic information was proposed (Refahi et al., 2016). Here we provide evidence that organ initiation in the SAM is indeed dependent on temporal integration of the auxin signal. Our quantitative analysis of the dynamics of auxin distribution and response supports a scenario in which rhythmic organ initiation at the SAM is driven by the combination of high-precision spatio-temporal graded distributions of auxin with the use of the duration of cell exposure to auxin, to temporally differentiate sites of organ initiation (Figure 6). Importantly our results suggest that a time integration mechanism is essential for rhythmic organ patterning in the SAM since auxin-based spatial information pre-specifies several sites of organ initiation and is thus unlikely to provide sufficient information (Video 1). Whether temporal integration of auxin information exists in other tissues remains to be established.

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 cell radial movement. Cells exiting the CZ that are exposed to high auxin levels progressively acquire competence for transcriptional response. This leads to activation of transcriptional responses with a delay close to the system period, the plastochron.

We provide evidence that temporal integration of the auxin signal likely requires the effectors of the auxin signaling pathway. Activation of transcription downstream of auxin by ARFs relies on chromatin remodeling, increasing the accessibility of ARF targets and possibly allowing for the recruitment of histone acetyltransferases (Wu et al., 2015), together with the release of histone deacetylases (HDACs) from target loci through degradation of Aux/IAA repressors (Long et al., 2006). Chromatin state change is one mechanism that allows the temporal integration of signals in eukaryotes, including plants (Angel et al., 2011; Coda et al., 2017; Nahmad and Lander, 2011; Sun et al., 2009). It is thus plausible that time integration of the auxin signal in cells leaving the CZ is set by progressive acetylation of histones triggered by ARFs at their target loci. As chromatin deacetylation also represses auxin signaling in the CZ (Ma et al., 2019), balancing the acetylation status of ARF target loci could provide a mechanism to tightly link stem cell maintenance to differentiation by precisely positioning organ initiation at the boundary of the stem cell niche, while at the same time allowing sequential organ initiation. Temporal integration might as well rely on mechanisms that fine-tune the intracellular distribution of auxin, such as auxin metabolism but also intracellular transport (Sauer and Kleine-Vehn, 2019). Determining how different mechanisms might act in parallel to provide a capacity to activate target genes as a function of auxin concentrations over time will require further analyses. It will notably be important to determine whether other ARF than ARF3 act in the temporal integration of auxin.

The existence of high definition spatio-temporal auxin gradients suggests that as for several morphogens in animals (Nahmad and Lander, 2011; Dessaud et al., 2007; Scherz et al., 2007; Maden, 2002) the robustness of SAM patterning results from highly reproducible spatio-temporal positional information. Our results indicate that auxin maxima could first emerge from the CZ at the confluence of centripetal auxin fluxes. Confluences creating auxin maxima would at the same time divert fluxes away from areas where auxin minima appear (Figure 3M). Our analysis raises the question of how auxin transport could generate this high definition signal distribution and whether the different models that have been proposed can explain this distribution (Bainbridge et al., 2008; Stoma et al., 2008; Jonsson et al., 2006; Smith et al., 2006a; van Berkel et al., 2013). Further analysis of the spatio-temporal control of auxin distribution needs also to consider that early developing flowers act as auxin production centers. These flowers could not only provide a memory of the developmental pattern through lateral inhibition but also contribute positively to a self-sustained auxin distribution pattern by providing auxin to the system (Figure 3M). Finally, our work indicates that the stem cell niche could act as a system-wide organizer of auxin transport, consistent with previous work (de Reuille et al., 2006). This could provide another layer of regulation tightly coordinating differentiation with the presence of a largely auxin-insensitive stem cell niche (Vernoux et al., 2011; Ma et al., 2019).

Materials and methods

Key resources table
Reagent type
(species) or resource
DesignationSource or referenceIdentifiersAdditional
information
Genetic reagent Arabidopsis thalianapPIN1:PIN1-GFP (Col-0)Benková et al., 2003
Genetic reagent Arabidopsis thalianapCLV3:mCherry-NLS (Col-0)Pfeiffer et al., 2016
Genetic reagent Arabidopsis thalianapYUC1-11:GFP
(Col-0)
Liu et al., 2016; Robert et al., 2013
Genetic reagent (Arabidopsis thaliana) yuc1 yuc4/+ pDR5rev::GFP (Col-0)Robert et al., 2013
Genetic reagent (Arabidopsis thaliana)ett-22 (Col-0)Pekker et al., 2005
Genetic reagent (Arabidopsis thaliana)pid-14 (Col-0)Huang et al., 2010
Genetic reagent Arabidopsis thalianapRPS5a:DII-VENUS-N7-p2A-TagBFP-SV40 (Col-0)This studyqDIIRequest to teva.vernoux@ens-lyon.fr
Genetic reagent Arabidopsis thalianapDR5rev:2x-mTurquoise2-SV40 (Col-0)This studyRequest to teva.vernoux@ens-lyon.fr
Chemical compound, drugTrichostatin AInvivogenmet-tsa-10.005 mM
Chemical compound, drugIndole-3-acetic acid sodium saltSigma-AldrichI51480.2, 1.0, 5.0 mM
Software, algorithmRStudioRStudio Team, 2015RRID:SCR_000432
Software, algorithmImage JRRID:SCR_003070https://imagej.net
Software, algorithmNumPyRRID:SCR_008633http://www.numpy.org
Software, algorithmSciPyVirtanen et al., 2020RRID:SCR_008058http://www.scipy.org
Software, algorithmVTKSchroeder et al., 2006RRID:SCR_015013http://www.vtk.org
Software, algorithmscikit-imagevan der Walt et al., 2014http://scikit-image.org
Software, algorithmsam_spaghettiThis study
Cerutti et al., 2020
https://gitlab.inria.fr/mosaic/publications/sam_spaghetti/
OtherPropidium iodide solutionSigma-AldrichP48640.1 mM

Plant material and growth conditions

Request a detailed protocol

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

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

Microscopy

Request a detailed protocol

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

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

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

Quantification and statistical analysis

Request a detailed protocol

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

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

Data and software availability

Request a detailed protocol

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

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

Appendix 1

Definition of a conceptual frame for models and analysis

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 (phyllotactic dynamical system)

Let a phyllotactic dynamical system S be a finite set of primordia considered over a time interval T=[tmin,tmax] and such that:

  • At every time t𝒯, each primordium p𝒮 is characterized by its current state {τp(t),xp(t),yp(t)} where:

    • τp[τmin,τmax] is called the developmental state of the primordium.

    • xp=[rp(t),θp(t),zp(t)]3 is the spatial position of the primordium in a cylindrical coordinate system, the origin of which is called the center of the system.

    • yp(t)d is a vector describing the physiological state of the primordium.

  • The developmental state τp is a continuous strictly increasing function of time. Note that consequently, for every p𝒮, τp is a bijection between 𝒯 and τp(𝒯).

  • In the case where 0τp(𝒯), the time t0,p=defτp-1(0) is called the initiation time of the primordium p.

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

    • (1) {X:[τmin,τmax]R3p𝒮,XpR3t𝒯,xp(t)=Xp+X(τp(t))Y:[τmin,τmax]Rdp𝒮,t𝒯,yp(t)=Y(τp(t))
  • 𝒮 is equipped with a strict total order denoted < that verifies:

    • (2) p,q𝒮,p<qt𝒯,τp(t)<τq(t)

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

Property 1

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

(3) r𝒮p<r<qdefq=p+1

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

Definition 2 (plastochron)

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

(4) p,p+1𝒮,t0,p+1=t0,p-T

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

Definition 3 (steady development)

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

(5) p𝒮,t𝒯,τp(t)=1T(t-t0,p)

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

Property 2

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

(6) p𝒮,t𝒯,τp(t+T)=τp(t)+1

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

Property 3

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

(7) p,p+1𝒮,t𝒯,τp+1(t)=τp(t)+1

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

Property 4 (developmental stage)

In a system S with a steady development of plastochron T, at any time tT, the rounding function of the developmental state:

(8) P(,t):𝒮[[Pmin(t),Pmax(t)]]ZpP(p,t)=τp(t)+12

is an isomorphism. We call (p,t) the developmental stage of primordium p at time t. If (p,t)=k, we say that the primordium p has the label k at time t.

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

Property 5

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

(9) p,p+1𝒮,t𝒯,(p+1,t)=(p,t)+1

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

Definition 4

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

(10) p,q𝒮,kt𝒯,(q,t)=(p,t)+kdefq=p+k

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

Definition 5 (regular phyllotaxis)

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

(11) p,p+1𝒮,Xp+1=Rz(α)Xp=Xp+(0α0)

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

Property 6 (spatio-temporal periodicity)

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

(12) p,p+1𝒮,t𝒯,{xp+1(t)=Rz(α)xp(t+T)yp+1(t)=yp(t+T)

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

Definition 6 (clear regular phyllotaxis)

In a system S with a steady development of plastochron T and a regular phyllotaxis of divergence angle α, if α is not a simple fraction of 2π, that isif:

(13) a,b[[1,|𝒮|]]α=ab2π

we say that S has a clear regular phyllotaxis.

Property 7 (ordering on a clear regular phyllotaxis)

In a system S with a steady development of plastochron T and a clear regular phyllotaxis of divergence angle α, then the primordia angles {θppS} are sufficient to determine the strict total order on S:

(14) p,q𝒮,t𝒯,kZθq(t)θp(t)=kα[2π]q=p+k

With this conceptual frame in mind, we can define thoroughly the general problem addressed when labeling the primordia on a meristem observation, that is typically when one tries to estimate where is 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.

Problem 1 (assignation of developmental stages)

Given a system S, observed at Kt discrete time points {tii[[0,Kt[[} in which, for every pS and for every ti, only a partially observed state {x~p(ti),y~p(ti)} is available;

Find for every tii[[0,Kt[[ an estimated developmental stage function ^(,ti) that verifies:

(15) p,p+1𝒮,^(p+1,ti)=^(p,ti)+1

and that minimizes the average staging error of the primordia:

(16) ϵP1|𝒮|p𝒮ϵP(p)whereϵP(p)=1Kti=0Kt1|P^(p,ti)τp(ti)|

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

Definition 7 (k-maintaining system)

We say that a system S is k-maintaining (k) if at all times, there is a primordium that has the label k:

(17) t𝒯,p𝒮,(p,t)=k

Property 8 (reduced k assignation problem in a clear regular phyllotaxis)

Let S be a k-mantaining system with a steady development of plastochron T and a clear regular phyllotaxis of divergence angle α. The solution to the assignation of developmental stages problem can be reduced to:

Find for every tii[[0,Kt[[, the primordium pS such that ^(p,ti)=k.

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

Definition 8 (k-characteristic state function)

Let S be a system with a steady development of plastochron T. We say that the physiological state function Y is k-characteristic (k) if there exists a value set Γkd such that:

(18) τ[τmin,τmax],Y(τ)Γk|τ-k|<12

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

Property 9 (k-characterization solution)

Let S be a k-mantaining system with a steady development of plastochron T and a clear regular phyllotaxis of divergence angle α and a k-characteristic state function. The reduced solution to the assignation of developmental stages problem given by:

(19) i[[0,Kt[[,^(p,ti)=ky~p(ti)Γk

has a staging error ϵ(p) bounded by 1/2 for every pS.

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

Property 10 (stage stationarity condition)

Let S be a system with a steady development of plastochron T. We say that a developmental stage assignation ^ is stationary if it is the same at all times of observation, that is if:

(20) ^*:𝒮p𝒮,t𝒯,^(p,t)=^*(p)

If S is observed during a time interval T=[tmin,tmax] smaller than its plastochron T, then there exists a stationary developmental stage assignation with a staging error bounded by 1/2 for every pS:

(21) tmaxtmin<TP^:𝒮Zp𝒮,ϵP(p)=1tmaxtmintmintmax|P^(p)τp(t)|dt<12

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

Property 11 (stationary k-characterization solution)

Let S be a k-mantaining system with a steady development of plastochron T and a clear regular phyllotaxis of divergence angle α and a k-characteristic state function. If S is observed during a time interval T=[tmin,tmax] smaller than its plastochron T, then the reduced stationary solution to the assignation of developmental stages problem given by:

(22) ^*(p)=k1Kti=0Kt-1y~p(ti)Γk

has a staging error ϵ(p) bounded by one for every pS:

(23) tKt1t0<Tp𝒮,ϵP(p)<1

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

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

Appendix 2

Effects of variability on a theoretical 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 plastochrons and divergence angles are drawn from random distributions centered on a common average value. By varying the levels of noise on both angular positions and plastochrons, we assess how variability impacts the overlapping of phyllotactic patterns at the scale of a population.

Let us consider a 2D phyllotactic dynamical system 𝒮 (see Appendix 1) formed by consecutive organ primordia observed on a temporal interval 𝒯. At every time t𝒯, each primordium labelled p,p[[0,pmax]] is represented by its developmental stage τp and by its coordinates in a 2D cylindrical reference frame:

(24) 𝒮(t)={(τp(t),rp(t),θp(t))p[[0,pmax]]}

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

(25) t𝒯,p]]0,pmax]],τp(t)=τp-1(t)+1

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

(26) t𝒯,p[[0,pmax]],{θp(t)=θ0+pαrp(t)=Reβτp(t)

This can be translated into incremental equations to obtain a recursive definition of the state of the system, knowing the state of the primordium labelled 0 at each time t𝒯:

(27) t𝒯,p]]0,pmax]],{θp(t)=θp1(t)+αrp(t)=rp1(t)eβ

We set ourselves in a context where all the considered phyllotactic systems have previously been aligned on 0, that is where t𝒯,θ0(t)=0 and τ0(t)[-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 plastochron time T

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

(28) t0,p-1-t0,pT𝒩(1,στ)

which translates into:

(29) t𝒯,τp(t)-τp-1(t)𝒩(1,στ)

Consequently, we can formulate the recursive definition of the system as the drawing of 2pmax random variables:

(30) t𝒯,p]]0,pmax]],{θp(t)θp1(t)𝒩(α,σα)rp(t)rp1(t)eβ𝒩(1,στ)

We simulate a population of such systems by generating K𝒮 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 (r0,θ0)=(eβτ0,0) and construct the system recursively by drawing the corresponding random variables. This way, we obtain a population of organ primordia positions identified by their rank p (Appendix 2—figure 1A).

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

(31) 𝒫={((ri,θi),pi)i[[0,(pmax+1)K𝒮[[}

To answer the separability question, we measure to which extent the identically labeled points 𝒫p={(ri,θi),pi=p} are separable by applying an unsupervised clustering algorithm. For this, we clustered the points using a k-means algorithm. The resulting clusters can be separated by linear boundaries in the Voronoi diagram associated with their centroids {(θp^,rp^)p[[0,pmax]]}. We measure the linear separability of primordia by looking if points with the same label gather inte the same Voronoi cell.

We use prior knowledge by setting the number of components of the k-means algorithm to pmax+1 and by setting the initial centroids {(θp0^,rp0^)p[[0,pmax]]} on the positions of the primordia in a model without any noise:

(32) p[[0,pmax]],{θp0^=pαrp0^=Reβp

After convergence, the algorithm returns pmax+1 centroid points that we use to construct the Voronoi diagram. This actually defines a predictor for the estimated primordium rank pi^ by looking inside which cell of the diagram lies a given point (Appendix 2—figure 1B).

(33) i[[0,(pmax+1)K𝒮[[,pi^=argminp[[0,pmax]]riθiri^θi^

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

(34) v=1(pmax+1)K𝒮i=0(pmax+1)K𝒮-1δpi^pi

If v equals to 1, it means that the primordium points group into perfectly identifiable clusters. This can be interpreted as the fact that K𝒮 randomly sampled individuals can be superimposed perfectly, once they have been centered and aligned on their primordium which is the closest to the 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 and 2 plastochrons, first with the values R=30μm, α=137.51 and β=0.23 corresponding to actual measured data (Figure 1—figure supplement 2D–E). As expected, increasing the angular variability creates more elongated clusters (Appendix 2—figure 1C) that still appear separated. Yet, the Voronoi separation introduces confusion between neighboring primordia, more specifically making p and p+3 overlap (Appendix 2—figure 1F). Interestingly, when we increase the plastochron variability (Appendix 2—figure 1D), the confusion concerns rather p and p+5 (Appendix 2—figure 1G). In both illustrated cases, the separability score drops markedly below 95%, while the separability of the actual observed data has been evaluated in the same way at 100% (Appendix 2—figure 1E).

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

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

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

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

Another interesting feature evidenced by this analysis is the influence of the motion speed of primordia separability. We explored the β×στ parameter space by fixing the value of σα to the observed one, and varying the motion coefficient β between 0 and 0.4 (Appendix 2—figure 1K). 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 plastochron 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.

Appendix 2—figure 1
A geometrical model of primordia distribution enables estimating the plastochron variability of the SAM.

(A) Primordium points are generated from a computational phyllotactic model with a control over the variability in angular positioning and plastochron duration. (B) Linearly separable clusters in the resulting point cloud are identified using an unsupervised algorithm with prior information. The obtained labeling in primordia ranks is compared with the theoretical one to compute a separability measure. (C–D) Increasing the azimuthal variability creates clusters that are more difficult to separate and generates confusion between p and p+3. (E) The primordia points from the observed experimental data form perfectly separable clusters (F–G) Increasing the plastochron variability creates clusters that are more difficult to separate and generates confusion between p and p+5. (H–J) Measured anglular and radial deviations of both observed and simulated primordia allow to determine plausible values for the angular varibility σα=3.6 (I) and plastochron variability σT=0.22 (J). The prmordia sdistriibution generated with this parameters (H) is the one that metches best the observed deviation of primordia clusters. (K) Separability evaluated by varying speed coefficient and plastochron variability. Modifying the radial speed of primordia changes the tolerance of the system to azimuthal and plastochron 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.

Appendix 3

Quantitative analysis of nuclei image signals

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

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

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

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

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

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

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

Appendix 3—figure 1
Automatic quantification pipeline for the time-lapse microscopy images of nuclei-targeted fluorescence signals.

To obtain quantitative data from the images produced under the microscope, various sequential processing steps need to be performed, from the extraction of the relevant objects (nuclei positions with their different channel intensity values) to the geometrical characterization and the spatio-temporal registration of the tissues, to finally get a complete, aligned and consistent dataset gathering all the imaged meristems.

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

Automatic nuclei detection

We consider that a 3D image consists of an array I of size Kx×Ky×Kz filled with values taken in an integer intensity interval . The elements of this 3D array are called voxels. In the case of a 16-bit-encoded unsigned integer image, the intensity interval =[[0,216-1[[. We denote:

(35) I={Iijk(i,j,k)[[0,Kx[[×[[0,Ky[[×[[0,Kz[[}.

The voxels in the array grid can be projected into a physical space Ω3 through a mapping function 𝐱:[[0,Kx[[×[[0,Ky[[×[[0,Kz[[Ω that associates the array indices with a discrete, evenly spaced, 3D lattice ΩIΩ. The voxel size v=(vx,vy,vz)3 defines a spacing of the lattice that is potentially different on each dimension:

(36) ΩI={𝐱(i,j,k)=(ivx,jvy,kvz)(i,j,k)[[0,Kx[[×[[0,Ky[[×[[0,Kz[[}.

This mapping allow to define the image as a function I:ΩI that associates a 3D point 𝐱=(x,y,z) in the image definition space ΩI with a fluorescence intensity value I(𝐱), such that (i,j,k)[[0,Kx[[×[[0,Ky[[×[[0,Kz[[,

(37) I(𝐱(i,j,k))=Iijk.

In the case of a multichannel image, we denote IS the image channel corresponding to the signal S. Nuclei are detected for each meristem acquisition independently using the pRPS5a:TagBFP (Tag) channel, which we denote ITag.

As a first approximation, cell nuclei can be considered to appear as a roughly spherical blob of fluorescence intensity in the image ITag, with a limitedly variable radius. To detect them we convolve the image with a sequence of Kσ 3D isotropic Gaussian kernels of increasing standard deviations Ωσ={σll[[0,Kσ[[}, with σ0<<σKσ-1. The response to this filtering is expected to be maximal for homogeneous spheres of intensity with a radius close to the standard deviation of the Gaussian kernel.

The sequence of filtered images can be seen as a 4D Gaussian scale-space (Lindeberg, 1994) resulting in a 4D image, which by extension we denote ITag:ΩI×ΩσR4. In this image, the fourth dimension σ corresponds to scale. If we note G(σ) the discrete Gaussian kernel of standard deviation σ, the scale-space transform we use is defined by, (x,σ)ΩI×Ωσ,

(38) ITag(𝐱,σ)=σ(ITag*G(σ))(𝐱).

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

A 4D point (𝐱,σ) is then considered a local maximum if its response ITag(𝐱,σ) is higher than a threshold Imin, and higher that all the neighbouring responses at all scales. More formally, if we note (𝐱,σ)={𝐱ΩI𝐱-𝐱<σ} the discrete ball of radius σ in the image lattice ΩI, the 4D point (𝐱,σ) is a local maximum if it realizes the maximal value of response over (x,σ)×Ωσ, in other terms:

(39) (x,σ) is a local\ maximum of ITag(x,σ)argmax(x,σ)(x,σ)×Ωσ(ITag(x,σ)).

Let us denote 𝒫 the set of points 𝐱ΩI corresponding to a local maximum of ITag whose response is strictly greater than Imin. Each point in 𝒫 corresponds to a detected nucleus, identified by an integer index n and associated with a single spatial position Pn in physical coordinates, so that we can write

(40) 𝒫={Pnn𝒩=[[1,|𝒫|[[}.

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

Nuclei signal quantification

Each image channel is quantified at the level of every detected nucleus. The signal intensity value is obtained by computing a weighted average of the channel intensity IS around the position of the nucleus. The signal images showing some local subcellular noise, the raw voxel value IS(Pn) might not be fully representative of the whole nucleus. We chose to use a distance-based Gaussian weight, of constant radius σ𝒩 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μm. Practically for the signal S, we first filter the image channel IS by a Gaussian kernel of radius σ𝒩 and retrieve the values at the voxel positions of all detected nuclei, so that n𝒩,

(41) Sn=(IS*G(σ𝒩))(Pn).

For example, the local level of expression of the CLV3 gene, imaged using pCLV3:mCHERRY in the channel ICLV3 would we quantified as CLV3n=(ICLV3*G(σ𝒩))(Pn). In the case of the ratiometric auxin sensor qDII, we combine the information of two fluorescence channels ITag and IDII to compute the ratio of estimated signals for each nuclei point, so that n𝒩,

(42) qDIIn=DIInTagn.

Extraction of L1 cells

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 based on the pRPSa:TagBFP channel as a 3D triangle mesh ={𝒱,𝒯}, where:

  • 𝒱 is a set of vertices

  • Each vertex v𝒱 is associated with a 3D position Mv3

  • 𝒯 a set of triangles defined by triplets of indices (v1,v2,v3)𝒱3

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

Computation of the surface mesh

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

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

Detection of L1 nuclei

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

(43) n(v)=argminn𝒩(Pn-Mv).

We define 1 as the subset of 𝒩 formed by indices n of nuclei points Pn that are the closest point to at least one vertex of .

Rigid time registration

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

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

Registered nuclei points

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

(44) Pn0=(j=i-10Rtjtj+1)Pn.

Non-linear time registration

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

We compute this new transformation as a dense vector field that maps two consecutive pRPS5a:TagBFP images, that have previously been rigidly registered into the coordinate system of I(t0). We denote {I(ti)Tag0i[[0,Kt[[} the sequence of such registered images.

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

The resulting vector field is actually a vectorial image 𝐔titi+10:ΩI(t0)3 defined over the same grid as I(t0). It contains at each voxel position a 3D vector measuring the local total deformation of the tissue to go from the current frame to the next one.

We also perform the backwards non-linear registration by computing another set of non-linear transformations between I(ti+1)Tag0 and I(ti)Tag0 in the exact same manner and obtain the vector fields 𝐔titi+10 mapping the next frame to the current one.

Nuclei cellular motion

The cellular motion between two consecutive sequence frames taken at times ti and ti+1 in the forward direction can be estimated using the registered nuclei points P(ti)0={Pn0n𝒩(ti)} and the transformation 𝐔titi+10 that maps the current frame into the next one by a vector field. Each nucleus n can then be assigned a local displacement vector by looking at 𝐔titi+10(Pn0).

We deduce speed vectors 𝐯n0 measuring the local speed of cellular motion by dividing the displacement vector by the time interval: n1(ti),

(45) 𝐯n0=1ti+1-ti𝐔titi+10(Pn0)

2D maps of epidermal signal

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

(46) η(r)=12-12tanh(k(r-R))

The continuous signal map S^ is then defined based on a point cloud 𝒫={Pn3n𝒩} and the associated signal values {Snn𝒩} as : 𝐱3,

(47) S^(𝐱)=1n𝒩η(𝐱-Pn)n𝒩η(𝐱-Pn)Sn

Note that, for any point 𝐱 where the total density H(𝐱)=n𝒩η(𝐱-Pn) equals , the signal map is not defined. To make the map outlines closer to their actual support, we consider that S^(𝐱) is defined only if H(𝐱)12. This constraint is equivalent to consider that the implicit surface obtained with the point cloud 𝒫 as generator and η as potential forms the boundary of the definition domain of the function S^.

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

Appendix 3—figure 2
2D continuous maps of epidermal signals.

(A) To build a continuous 2D map, we diffuse the signal in space by computing a local average of discrete signal values using a kernel function whose extent and sharpness are set by two parameters R and k. (B). Optimal values of these parameters were determined on the signal of main interest qDII and 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.

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

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

Sequence SAM reference frame determination

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

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

  • the position 𝐜=(x𝐜,y𝐜,z𝐜)3 of the apex center in the central zone (CZ) of the meristematic dome

  • the unitary vector 𝐚3 of the main radial symmetry axis of the meristematic dome

  • the unitary radial vector 𝐫3 ( 𝐚𝐫) of the direction of the last initiated organ primordium (labelled 0) relatively to 𝐜

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

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

Detection of the CZ center 𝐜

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

In the resulting 2D map CLV3^, the CZ appears as a wide isotropic peak of signal intensity (Appendix 3—figure 3D). We assume that the CZ is the largest area of high CLV3 intensity in the tissue we consider. We extract a central zone domain by thresholding the 2D map CLV3^ and keeping the largest connected component, from which we compute the 2D center (x𝐜,y𝐜) and the area A𝐜.

To make the estimation more robust, we perform this CZ extraction using a range of K𝐜 threshold values that depend on the average signal intensity of the sequence {ckCLV3¯k[[0,K𝐜[[}. In the end, we estimate (x𝐜,y𝐜) as the average 2D center and A𝐜 and the average area of all these central zone domains obtained with K𝐜=7 and ck=1.2+0.1k.

The resulting distribution of estimated CZ radii r𝐜=A𝐜π among the considered individuals showed a very peaked distribution around an average value of 28 µm (Appendix 3—figure 3E).

Finally, the third coordinate z𝐜 of the 3D center of the CZ is computed as the estimated value z0^ of the z coordinate of L1 nuclei points in 0, taken at the center point (x𝐜,y𝐜) detected above.

Appendix 3—figure 3
Landmark-based alignment of SAMs using 2D continuous maps.

(A) Seen in 3D, the SAM shows a dome-like structure surrounded by a spiral of primordia. (B). The center and the main axis of the dome, the direction of the last intiated primordium and the orientation of the phyllotactic spiral are key landmarks of the SAM geometry. (C). Knowing their position allows to define a reference frame * in which all individuals can be superimposed. (D). The 2D projected map of epidermal CLV3 signal allows locating precisely the center 𝐜 and main axis 𝐚 of the meristem and estimating the extent of the CZ. (E). 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). (F). The 2D projected map of epidermal qDII signal allows detecting the direction 𝐫 of the 0 primordium as the maximal concentration of auxin in the PZ. The phyllotactic orientation o is set manually. (G) The rigid transformation is applied to the detected nuclei to align the individuals in the reference frame *. (H) The same is done with images and a curved L1 slice is used to display only epidermal cells.

Estimation of the vertical axis of the apex 𝐚

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

This vector is characterized by a rotation matrix R𝐚 to allows to transform the reference frame 0 of the image to a new reference frame 𝐚 where the origin is 𝐜 and the z axis corresponds to 𝐚.

We estimate the rotation matrix by minimizing the tilting of the surface of the CZ in this new reference frame. To do so, we look at the dispersion of z coordinates of the nuclei expressed in this reference frame as a function of their radial distance. We minimize the error between the z coordinate of the L1 nuclei points and their local average z¯. To take into account the extent of the meristematic dome, this error is weighted by a Gaussian function of the distance of points to 𝐜 of standard deviation σr=20μm.

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

Detection of the 0 direction 𝐫

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

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

Determination of the phyllotactic orientation

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

(48) Ro=(1000o0001).

Aligned L1 nuclei points

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

(49) Rt0=RoRrRaco1

We note 𝒫*={Pn*n1} the positions of first-layer nuclei points in this common reference frame, defined by, n1,Pn*=R*t0Pn0. From now on, we will always consider that the nuclei points have been aligned in the common SAM reference frame *, and that the signal maps are computed using the aligned positions 𝒫*.

Aligned L1 image slices

We also use the resulting alignment transformation to register the original images into the common reference frame by applying the same transform R*t0 to the registered image I(ti)S0. This creates a registered image I(ti)S* expressed over a new voxel grid ΩI* that is centered on 0 in x and y and slightly shifted in z so that:

(50) ΩI*=ΩI-(Kx2vxKy2vyKz8vz).

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

(51) I(ti)S*(x,y)=I(ti)S*(x,y,z*^(x,y)vzvz).

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 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 (both maxima and minima) in this 2D map and in labeling them in a second time by organ primordium rank.

Detection of extremal regions in the auxin map

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

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

(52) Λu={1if{uqDII^=0uuqDII^>00otherwiseΛu+={1if{uqDII^=0uuqDII^<00otherwise

If we consider all possible directions, defined by their angle ω[-π,π], a valley point (respectively a ridge point) is a point that achieves local minimality (respectively local maximality) along a large proportion of unitary vectors 𝐮(ω)=(cos(ω),sin(ω)). Consequently, we define the scalar fields Λ-:2[0,1] and Λ+:2[0,1] measuring respectively the global valleyness and ridgeness of a point in the field qDII^ as:

(53) Λ-=12π-ππΛ𝐮(ω)-𝑑ω Λ+=12π-ππΛ𝐮(ω)+𝑑ω.

We first look for saddle points in order to eliminate them and disconnect the otherwise continuous networks of ridges and valleys. Saddle regions 0 are formed by points that belong to both a valley and a ridge and are detected using the geometric mean of valleyness and ridgeness fields, then valley regions - and ridge regions + can be found outside the saddle regions, using different thresholds λ0, λ+ and λ- in each case:

(54) 0={xR2Λ+(x)Λ(x)λ0}={xR20Λ(x)λ}+={xR20Λ+(x)λ+}

We perform this region extraction with the values λ0=0.01 and λ+=λ-=0.03. Then we consider the union of the connected components of each set 0, - and +, which results in a set of connected regions 𝒞={cll[[0,K𝒞[[} as the ones delineated in Appendix 3—figure 4B. These regions are guaranteed to have no overlap as long as λ0<λ+ and λ0<λ-.

Description of extremal regions

Each of these connected extremal regions can potentially match the auxin accumulation or depletion zone characterisitic of a given primordium. In order to determine this matching we characterize the information contained in each of the extremal regions by defining a set of descriptor functions:

  • type:𝒞{-1,0,1}: whether the region is a connected component of -, 0 or + respectively.

  • pos:𝒞2: a single spatial position corresponding to the point of extremal value of the qDII^ field.

  • area:𝒞: the area of the connected region.

  • extr:𝒞[0,1]: the maximal value of extremality (valleyness, geometric mean or ridgeness).

  • qDII:𝒞: the extremal value of the qDII^ field in the connected region.

For each region, we define two significance scores γ-:𝒞[0,1] and γ+:𝒞[0,1] measuring how good a candidate the region would be respectively for a minimal or a maximal area of the qDII^ map. Each of these two scores is computed using the following elementary score functions, which only use the region descriptors:

  • A central zone exclusion score γCZ impeding the regions from lying within the central zone of the meristem, defined as:

    • (55) γCZ(cl)={1ifpos(cl)>r𝐜20otherwise.
  • An area score γarea that is closer to one when the considered region is large, defined as:

    • (56) γarea(cl)=1-1area(cl).
  • An extremality score γext that is closer to one when the considered region has a large value of extremality, defined as:

    • (57) γextr(cl)=min(12+extr(cl),1).
  • A signal value score, γqDII that is closer to one when the value of qDII^ is low for minimal regions and high for maximal regions (saddles being considered as maximal regions of lower signal), defined as:

    • (58) γqDII(cl)=(1-qDII(cl))δtype(cl),-1+(0.2+qDII(cl))δtype(cl),0+(qDII(cl))δtype(cl),1
    • where δ is the Kronecker symbol.

We compute the significance scores γ- and γ+ as a multiplicative combination of these elementary scores to obtain maximal values respectively for large, marked valleys of low qDII^ value and large, marked ridges or marked saddles of high qDII^ value outside the middle of the CZ (Appendix 3—figure 4C). Therefore we define the significance scores of extremal regions as:

(59) γ(cl)=γextr(cl)γCZ(cl)γqDII(cl)(δtype(cl),0+γarea(cl)δtype(cl),1).γ+(cl)=γextr(cl)γCZ(cl)γqDII(cl)γarea(cl)δtype(cl),1
Appendix 3—figure 4
Detection of primordia as extremal regions of auxin.

(A) The auxin distribution forms a landscape where accumulation froms troughs or valleys in the qDII map, whereas depletion forms peaks or ridges. (B). Extremal regions (valleys in purple, ridges in orange and saddles in green) are identified from the qDII map. (C). Each region is associated with a unique spatial position and characterized by a significance score. (D). A geometrical model of primordia distribution is used to compute a confidence score for each extremal region, relatively to each considered primordium rank, depending on its spatial postion. (E). The primordium rank k is assigned to the extremal regions with the highest combination of significance and spatial confidence score relatively to k (F). This allows to identify the auxin accumulation zoned and the potential auxin depletion zones associated with each primordium in the considered SAM. (G). 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.

Geometrical model of primordia organization

To help predict the position of primordia, we used a simple a priori geometrical model of primordia spatial organization in a 2D projected space. This model associates each point (r,θ)2 expressed in polar coordinates with a score γk(r,θ) reflecting the confidence of finding the kth primordium at the position (r,θ).

This score is decomposed into an angular score γk(θ) and a radial score γk′′(r) and such that γk(r,θ)=γk(θ)γk′′(r). The angular score makes the assumption that pimordia are organized in a regular spiral of divergence angle α*=2πϕ2137.5, and is therefore close to one when the angular coordinate of a 2D point is close to the theoretical angle kα* of the kth primordium, with an angular tolerance Δθ(k) (Appendix 3—figure 4D), so that:

(60) γk(θ)=cos(min(π|θ-kα*|Δθ(k),π2)).

The radial score makes the assumption that primordia move radially so that they are more likely to be found in a range of radial distances [rmin(k),rmax(k)]+ that gradually increases with primordium rank k (Appendix 3—figure 4D). We define the radial score γk′′, so that it is equal to one within the range and decreasing outside, by:

(61) γk′′(r)=min((rmax(k)r)2,1)min((rrmin(k))2,1).

We define the radial distance range as a relative measure on the CZ radius r𝐜, so that rmin(k)=ρmin(k)r𝐜 and rmax(k)=ρmax(k)r𝐜. The relative distances ρmin and ρmax are affine functions of the primordium rank i, such that:

(62) ρmin(k)=ρmin0+ρmin1kρmax(k)=ρmax0+ρmax1k.

The product of these two scores results in the delineation of smooth windows in the 2D space corresponding to the areas of high likelihood of encountering the primordium k (Appendix 3—figure 4D).

Assignment of primordia to extremal regions

We consider that primordia are characterized by the following specifications:

  • A primordium p is necessarily associated with the existence of a valley in 𝒞 and possibly with the additional existence of a ridge or a saddle in 𝒞.

  • A primordium p must respect several conditions:

    • The valley associated with p must locally be the most significant minimal region, that is have a locally maximal value of γ-.

    • If p is labelled k, the valley point associated with p must have a high spatial confidence score γk.

    • The ridge or saddle point associated with p must lie between the center of the meristem and the radial position of the valley point.

    • The ridge or saddle associated with p must locally be the most significant maximal region, that is have a locally maximal value of γ+.

    • If p is labelled k, the ridge or saddle point associated with p must have a high spatial confidence score γk.

For each primordium rank k[[kmin,kmax]], we make the decision of assigning the label k to a valley region cl𝒞, if the region achieves the highest combined score γk-(cl)=γk(pos(cl))γ-(cl) among the extremal regions, and that this score is greater than a threshold γmin:

(63) l-(k)={argmaxl[[0,K𝒞[[(γk-(cl))ifmaxl[[0,K𝒞[[(γk-(cl))>γminotherwise.

When the valley point exists (l-(k)) we use its position to restrict the search of a ridge or saddle point to the area located between the center the center and the valley region cl-(k) by defining an additional score function γk±:𝒞{0,1} as:

(64) γk±(cl)={1ifpos(cl)<pos(cl-(k))0otherwise.

In a second time, we make the decision of assigning the label k to a ridge or saddle region, if it achieves the highest combined score γk+(cl)=γk(pos(cl))γ-(cl)γk±(cl) among the extremal regions, and that this score is greater than Γmin:

(65) l+(k)={argmaxl[[0,K𝒞[[(γk+(cl))ifmaxl[[0,K𝒞[[(γk+(cl))>γminotherwise.

Each identified primordium extremal point is then associated with a spatial position (ri-,θi-) or (ri+,θi+) (Appendix 3—figure 4F) which can be used to compute estimates of all the spatial signals (starting with qDII^) and to track signals at the level of organ primordia.

Evaluation of primordia detection

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. The assignment of primordium ranks was made on the detected extremal regions 𝒞 so that manual and automatic assignment could be compared quantitatively.

This manual labeling allowed us to perform an evaluation of the detection method at the scale of the whole population of meristems (N = 63). We evaluated the performance by counting:

The number VP of correctly assigned primordia extremal points (corresponding to the case where the primordium rank k has been assigned to an extremal region in the automatic method whereas it has been assigned to the same one in the manual labelling)

  • The number MP of falsely assigned primordia extremal points (corresponding to the case where the primordium rank k has been assigned to an extremal region in the automatic method whereas it has been assigned to a different one in the manual labelling)

  • The number FP of falsely detected primordia extremal points (corresponding to the case where the primordium rank k has been assigned to an extremal region in the automatic method whereas it has not been assigned to any region in the manual labelling)

  • The number FN of falsely undetected primordia extremal points (corresponding to the case where the primordium rank k has not been assigned to any region in the automatic method whereas it has been assigned to an extremal region in the manual labelling)

The measure we used is the Jaccard index j=VPVP+MP+FP+FN and our evaluation lead to an overall average of 75.4% correct detection, rising to 87.3% when considering only primordia ranks between -2 and 3 (Appendix 3—figure 4G). More precisely, both the auxin accumulation and depletion zones of primordia 0, 2 and 3 are remarkably well detected (nearly 95% correct detection).

The loss of performance at 1 comes from the common missed detection of the depletion zone that only starts forming at this stage, and is only visible as a saddle point in the qDII^ map. For the rest, most of the errors are at early stages where auxin accumulation is not very marked or at later stages where the organs get larger, the auxin landscape more complicated (leading to a large number of false assignments) and curvature more important (making the 2D projection less relevant, especially at the tip of the organ where auxin accumulation takes place).

Data presented in Figure 1G, Figure 2D–E, Figure 3J, Figure 1—figure supplement 2D and Figure 1—figure supplement 3 were obtained using the manually corrected positions of organ primordia.

L1 dynamic signal maps

Developmental state estimation

We assume that all the observed meristems develop at a comparable rate of 2 new organs per day, corresponding to a plastochron of 12 hr, and that they can, at the first order, be considered synchronous (similarly labeled organ primordia being at comparable developmental stages. Consequently we defined a developmental state indexation as:

(66) τi=ti12,

where i[[0,Kt[[ and t0=0h.

Dynamic signal maps

At the scale of one time-lapse sequence, for which we have computed the 2D aligned signal maps {S^ii[[0,Kt[[} we approximate the continuous 2D+t signal using a density function of developmental state η defined by:

(67) η(δτ)=12-12tanh(kτ(δτ-Rτ)).

At any developmental state position τ, the estimated map is a weighted average of single time maps, where the weights are time-distance-based density coefficients computed using the ητ density function computed with a time radius Rτ and a time slope kτ:

(68) S^(r,θ,τ)=1i=0Kt-1ητ(|τi-τ|)i=0Kt-1ητ(|τi-τ|)S^i(r,θ)

Temporal extrapolation of dynamic maps

We evidenced that looking p plastochrons further in time is equivalent to rotating the system of the angle θp corresponding to the direction of the primordium p (see Figure 1—figure supplement 2G-H). We use this the spatio-temporal periodicity of the system to extrapolate the evolution of signals beyond the duration of acquisitions.

Considering a range of primordia stages [[pmin,pmax]], we use the angles θp- estimated with 2D primordium detection to apply rotations on the computed maps and derive a dynamic map that covers a larger temporal range as:

(69) S^(r,θ,τ)=i=0Kt-1p=pminpmaxητ(|(τi+p)-τ|)S^i(r,θ+θp-)i=0Kt-1p=pminpmaxητ(|(τi+p)-τ|).

Ultimately, if instead of using all time points of a given sequence, we use the maps of all time points from all the sequences in the previous formula, we reconstruct the average map S^ over a population of meristems. Such a map combines all the quantitative spatio-temporal information into one dynamic map reflecting the canonical behavior of the system. The result is presented for the Auxin signal in Video 1, where it was obtained with the parameter values Rτ=0.1, kτ=2, pmin=-3 and pmax=5.

Implementation details

  • Microscopy image preparation

    • The confocal images saved as CZI files through the ZEN software (ZEISS International, 2012) of the LSM-710 microscope are opened using a Python script (Gohlke, 2012) and split into independent channels that are saved separately as INR image files. This operation preserves all the information contained in the raw image format. In the specific case of acquisitions for which both pPIN1:PIN1-GFP (PIN1) and pRPS5a:DII-VENUS (DII) 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.

    • 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 six visible organs. This binary mask is then applied on all the channels, and masked channels are saved in separate INR files.

  • Nuclei detection and quantification

    • The nuclei detection algorithm has been implemented in the tissue_nukem_3d Python library using the Gaussian filtering functions provided by the SciPy library (Virtanen et al., 2020) for the construction of the scale space.

    • The surface mesh used for the extraction of L1 nuclei was computed using the implementations provided in the VTK (Visualization ToolKit) library (Schroeder et al., 2006) for the Marching Cubes algorithm, the mesh smoothing and the quadric decimation (vtkImageMarchingCubes, vtkWindowedSincPolyDataFilter, vtkQuadricDecimation). A version of the isotropic remeshing algorithm has been implemented in the cellcomplex library.

  • Image registration

    • The image registration operations are performed using the blockmatching computing library that comes embedded in the timagetk (Tissue Image ToolKit) image processing library dedicated to 3D microscopy images of multicellular tissues. Registered images are generally saved channel by channel in separate INR files, as well as the vector fields computed in the case of non-linear registration.

  • Continuous map generation

    • The methods for generating continuous maps of epidermal signal, and to detect extremal regions on such maps were implemented in the tissue_nukem_3d Python library using the Gaussian derivative functions provided by the SciPy library (Virtanen et al., 2020) for the computation of gradients.

Appendix 4

Quantitative analysis of PIN1 polarity

To study the cellular polarities of the auxin efflux carrier PIN1, we quantify information at the interface between cells, and the resulting data is best expressed as vectors. Therefore, it requires a whole different set of computational steps to analyze PIN1 polarities, and integrate polarity information at increasing levels:

  • Cell interface level: quantification of interface polarity from image intensities:

    • To quantify PIN1 polarity, we need to consider pairs of neighboring cells, and it is therefore necessary to reconstruct cell adjacency information from the images. The use of images with a cell wall staining allows partitioning the image grid into contiguous regions of voxels representing cells, which will give us neighborhood information as well as a way to reconstruct the geometry of cell interfaces. Based on these geometries, we need to compute a local PIN1 polarity at the scale of cell interfaces using the local signal intensity of the image.

  • Cell level: integration of cell interface polarities into cell polarity vectors:

    • To match the rest of the data that is defined at cell level, we need to convert the scalar information of polarization for each cell interface into a cellular information. The directional contributions from the different interfaces surrounding a cell necessarily result into a vectorial data at the level of the cell. Then, by applying the geometrical transform estimated as in Appendix B, we have to align not only the spatial positions of cellular objects but also the directions of cell polarity vectors.

  • Tissue level: computation of a continuous vector map of polarities that allows population averaging:

    • Finally, in order to identify global trends in the PIN1 polarity patterns across individuals, we need to use local averaging to build a continuous representation of signal that allows point-wise comparison, using the 2D map formalism. The only difference is that, in the case of polarities, the data we average consists of vectors. The result is therefore a continuous vector map on which we can compute other scalar properties, such as norm or divergence.

This alternative pipeline is illustrated in Appendix 4—figure 1, and it results in the continuous polarity maps shown in Figure 3I and Figure 3—figure supplement 2. Its development led to the introduction of an original method for the quantification of cell interface transporter polarity from images in 3D, which has been evaluated using super-resolution (using radial fluctuation [Gustafsson et al., 2016]) acquisitions of the same tissues (see Evaluation of cell interface polarity estimation). The different computational steps involved in the pipeline are detailed in the following sections:

Appendix 4—figure 1
Automatic quantification pipeline for the time-lapse microscopy images of PIN1 auxin efflux carrier.

The quantitative estimation of PIN polarity relies on the analysis of cell wall and membrane-marker images that need to be processed in a different way than nuclei marker images. An alternative automatic pipeline performs the necessary steps, from the segmentation of the cells and the extraction of L1 cell anticlinal interfaces to the quantification of signal distribution at each cell interface and the reconstruction of the polarized cell network.

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 cell wall segmentation step, we performed:

  • An adaptative histogram equalization for all z-slices of the membrane stacks to improve and normalize the contrast (Pizer et al., 1987);

  • 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 perform Gaussian smoothing and obtain smoother segmentation along the z axis (Fernandez et al., 2010);

  • Gaussian smoothing of the cell wall intensity image, with σ=0.2μm to reduce noise in the image when performing watershed segmentation (Fernandez et al., 2010);

  • Create a seed image from the nuclei coordinates to initialize the watershed algorithm (Fernandez et al., 2010);

  • Run the seeded-watershed algorithm with isometric smoothed intensity image and seed image (Fernandez et al., 2010).

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 Iseg that assigns an integer label to every voxel of the image grid ΩI on which the image IPI to segment is defined (see Appendix B). 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 Iseg). The background corresponds to a specific label, that is systematically set to one to ensure consistency between images. Each cell labeled c[[1,Kc]] is the then represented by a connected region Γc so that:

(70) Γc={𝐱ΩIIseg(𝐱)=c}

Cell barycenter extraction

To obtain the cells barycenter, we estimate the position Pc of the center of the cell labeled c in the segmented image Iseg as the average of the voxel coordinates of the cell region Γc, by c[[1,Kc]]{1},

(71) Pc=1|Γc|𝐱Γc𝐱.

Extraction of L1 segmented cells

Definition of cell adjacency

Adjacency relationship is defined through the notion of surface of contact between two regions in the segmented image. In a first time, we consider the 6-connectivity to define neighborhoods at voxel-scale. However, to be robust to potential segmentation errors, we estimate the area of all surfaces of contact representing walls between two cells c and c, and consider as neighbors only those for which the area is greater than a given threshold Amin.

We compute the area of contact A(Γc,Γc) as the sum of areas of the contact rectangles between pairs of 6-adjacent voxels so that one belongs to Γc and the other to Γc. An analysis of the distribution of wall areas (data not shown) led us to consider that a value of Amin=5μm2 would be suitable. We note 𝒜c the cells that are adjacent to the cell c, namely the labels that verify the surface of contact condition:

(72) 𝒜c={c[[1,Kc]]{c}A(Γc,Γc)>Amin}

Estimation of L1 layer

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. It is the definition of the epidermis to be in contact with the outside world, and all the regions of cells belonging to the 1 set should therefore be adjacent to Γ1. In some problematic cases, the segmentation algorithm can produce a background region that reaches the inside of the tissue, but we assumed that it was not the case in our images and that we could consider that 1=𝒜1.

Meshing L1 anticlinal cell interfaces

In order to quantify the PIN signal intensity at the level of each individual cell-cell surface of contact, we need to have a precise identification of the cell interfaces 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 Iseg, it will generally be too sensitive to noise and to image resolution to be used as such. We chose therefore to use a triangular mesh representation with a high resolution to represent accurately the cell interfaces.

To obtain such meshes, we apply the Marching Cubes algorithm (Lorensen and Cline, 1987) to each cell labelled c and represented by its connected region Γc of identically labeled voxels in the segmented membrane image Iseg. This produces a triangular mesh c=(𝒱c,𝒯c) where:

  • 𝒱c is a set of vertices

  • Every vertex v𝒱c is associated with a 3D coordinate Mv3

  • 𝒯c is a set of triangular faces t=(v1,v2,v3)𝒱c3 linking those vertices.

The cell meshes are generally closed (except on image borders) and have a voxel-like resolution.

Using Marching Cubes ensures us that, in an 8-voxel cube with only two labels c and c (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 c and c. We use this property to construct the cell interface mesh c,c=(𝒱c,c,𝒯c,c) of the interface between Γc and Γc as the restriction of one of the two meshes to the set of common vertex points:

(73) {𝒱c,c={v𝒱cv𝒱c,Mv-Mv<ϵ}𝒯c,c={t𝒯cvt,v𝒱c,c}

Each interface mesh undergoes then a phase of triangle decimation (Garland and Heckbert, 1997b) and isotropic remeshing (Botsch and Kobbelt, 2004b) 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 𝐧cc(v) at each vertex v of c,c, ensuring they all point from c to c, and the area of each triangle that allows us to estimate the total area Ac,c=Ac,c of the interface between cells c and c. Note that the Marching Cubes intersection, along with the decimation and smoothing, will produce a mesh that is smaller that the actual cell interface (the intersecting part does not extend to cell edges) and the interface area will be underestimated. On the other hand, the voxel-based area estimation A(Γc,Γc) is known to be largely overestimating, which ultimately provides a way to have both lower and upper bound estimates of the interface areas.

Quantification of PIN polarity at cell interfaces

PIN proteins are polarly distributed in cells. This potentially induces that on both sides of a given interface between two cells c and c, different concentrations on PIN can be observed. This suggests that a flow of auxin could be oriented from c to c if more PIN transporters are located at this interface in the cell membrane of c than in c.

This cell interface polarity orientation is denoted pcc=-pcc[-1,1] and is equal to 1 (respectively −1) when there exists a marked positive (respectively negative) difference between the PIN membrane concentrations of c and c at their interface. If the case where there is no difference, pcc=0, and a value between 0 and 1 reflects intermediate difference levels.

Appendix 4—figure 2
Estimation of PIN polarity at cell interface level.

(A) The triangular mesh representing the cell interface is used to generate a set of 3D cylinders locally orthogonal to the interface and placed at each vertex of the inside of the interface. (B) The radius of the cylinders is close to the typical resolution of the mesh. (C) For each cylinder, all the neighboring image voxels are projected onto its main axis, and kept if within the distance defined by the cylinder radius. Distance of the projected voxel to the interface mesh vertex is used as abscissa for the evaluation of 1d polarity (D) The set of all interface cylinders allow sampling the PIN image signal on either side of the cell-wall at different locations. (E) The 1-dimensional distributions of both PI and PIN image signals along each cylinder are used to locate precisely the cell-wall position and quantify signal levels on either side. (F) PIN levels are estimated left and right of the detected cell-wall abscissa up to a fixed distance of 0.6µ m. (G). Significant difference between left and right distributions of PIN levels across the interface allows deciding for local PIN polarity at the scale of the cell interface.

To estimate this polarity, we aim to sample the PIN signal on each side of the cell-wall separating c and c. For this, we proceed in several steps:

  • On each vertex v of c,c, we position a cylinder of radius r𝒞 and height 2d𝒞, centered on Mv and whose main axis is defined by the local normal vector 𝐧=𝐧cc(v) (Appenix 4—figure 2B-D). We excluded vertices located on the contour of the mesh where the estimation of 𝐧 is less robust.

  • With each voxel 𝐱ΩI in the cylinder, we associate a signed abscissa zv(𝐱)=𝐱-Mv,𝐧 (Appenix 4—figure 2C) that allows positioning it on the 1D axis of the cylinder.

  • Due to processing errors, the actual cell-wall (marked by PI) may be shifted from the position zv=0 where the cylinder intersects the interface mesh. We then consider the spatial distribution of PI intensities along the cylinder axis and model it using a parametric Gaussian function:

    • (74) PI(zv)=PI0+PI1e-(zv-zPI)2σPI2.
    • We estimate the parameters of this function by a least-squares optimization algorithm and use the value of zPI as a reference defining the actual position of the cell-wall (Appenix 4—figure 2E).

  • Finally, we compute the average of PIN intensities of voxels 𝐱 such that zPI-zmax<zv(𝐱)<zPI (respectively zPI<zv(𝐱)<zPI+zmax) to estimate the value PINcc(v) (respectively PINcc(v)) of PIN signal in cell c (respectively c) at the level of vertex v. Note that the definition is symmetrical, meaning that the result would be exactly the same if c and c were permuted.

By performing this two-sided estimation on every cylinder over the triangular mesh representing the cell interface, we end up with two paired signal distributions {PINcc(v)} and {PINcc(v)}. We test statistically if one side of the interface bears significantly more signal that the other by performing two one-sided Wilcoxon tests, a non-parametric version of the T-test for paired samples. They assess respectively if the means of the distributions (noted PINcc and PINcc) can be ranked, that is if PINcc>PINcc and if PINcc>PINcc respectively. The tests provide two p-values pvalcc and pvalcc based on which we make a decision on the existence of a polarity. We do not use a binary polarity value but define a polarity value pcc whose sign corresponds to the direction of polarity and absolute value lies among 1, 0.5, 0.25 and to account for situations where some uncertainty remains from the statistical tests:

(75) pcc={{1ifpvalcc<0.10.5ifpvalcc<0.15if(pvalcc<0.25)and(pvalcc>0.25)0.25if0.15pvalcc<0.25{1ifpvalcc<0.10.5if0.1pvalcc<0.15if(pvalcc<0.25)and(pvalcc>0.25)0.25if0.15pvalcc<0.250otherwise.

Note that the definition is symmetrical ensuring that pcc=-pcc. We also define the intensity of PIN signal at the level of the cell interface PINc,c=PINc,c as the average of the signal medians of either side, the barycenter Pc,c of the interface mesh, and interface normal vector 𝐧cc=-𝐧cc computed as the average of the non-contour interface vertex normals.

Evaluation of cell interface polarity estimation

We assessed whether this automatic polarity estimation, performed on confocal microscopy images with a typical resolution of 0.2μm (far beyond cell-wall thickness), manages to retrieve a difference of carrier concentration at the level of plasma membranes. To do so, we used super-resolution imaging to provide a set of ground truth polarities that we use to evaluate the result of our automatic method.

Super-resolution images were obtained using the Super-Resolution Radial Fluctuations (SRRF) technique (Gustafsson et al., 2016), which produces a 2D image with a resolution of 0.035μm (Appendix 4—figure 3A). The same tissue was immediately (0.5 to 2 hr) imaged under the confocal microscope to make sure that the local polarities will be comparable. The confocal image stack was then processed through the automatic pipeline described above. This provides notably unique identifiers to cells and allows to identify the visible cell interfaces by pairs of cell labels that can be matched to the outcome of the pipeline.

Individual cell interfaces were manually annotated by two independent experts in the ImageJ software using the same protocol. We designed a procedure to determine carrier polarities on 2D slices co-imaging membrane-located fluorescent carrier proteins and cell-wall staining inspired from Shi et al. (2017). The following steps are repeated for each cell interface, identified by a pair of cell labels coming from the automatically segmented PI image:

  • Draw a minimal amount of 3 lines using the ‘Straight line’ tool:

    • Going always from the same cell (left) towards the other (right)

    • Keeping them orthogonal to the apparent interface plane (Appendix 4—figure 3B)

  • For each line, sample the PI signal and the PIN1 signal along the line using the ‘Plot profile’ tool

    • With a line width set to 5 to 20 pixels for local averaging

  • Determine if the peak of PIN1 signal lies clearly and consistently on one side of the peak of PI signal.

    • If PIN1 if clearly and consistently to the left of PI, assign the polarity 1.

    • If PIN1 if consistently to the left of PI, but not clearly, assign the polarity 0.5.

    • If PIN1 if clearly and consistently to the right of PI, assign the polarity −1.

    • If PIN1 if consistently to the right of PI, but not clearly, assign the polarity −0.5.

    • If PIN1 is not consistently on one side of PI, or not clearly, assign the polarity .

Repeating this manual assessment for every pair of neighbor L1 cells visible in the SRRF image plane, each expert filled a data sheet using the same cell labels as in the segmented PI image and the manually determined polarity values. Among the expertized interfaces, we retained only those that were assessed by both experts and that were found in consensus, that is where either:

  • Both experts assigned the same polarity value

  • One expert assigned the value 1 (resp. −1) and the other 0.5 (resp. −0.5).

In the latter case, the ground truth polarity value was corrected to 1 (resp. −1), acknowledging the certainty given by one expert when both agree on the direction of polarity. This resulted in a set of consensus interface polarities as displayed in Appendix 4—figure 3C. In total, these consensus interfaces made up for more than 85% of the interfaces analyzed by both experts as shown in Appendix 4—figure 3D, and less than 4% of the interfaces found the experts in disagreement.

We compared these ground truth polarities with the polarity values pcc predicted by the automatic method for the same pairs of cells. For the sake of evaluation, we considered the two lower levels of certainty in the statistical test (corresponding to values 0.5 and 0.25) as the same uncertain value 0.5 in the ground truth assessment. We considered the cell pairs of the expertized interface to be in the order of their assigned polarity, that is such that all the ground truth polarity values p* are either , 0.5 or 1, which is possible due to the symmetric nature of the polarity values. The comparison translates into a table {C(p*,p)p*{0,0.5,1},p{-1,-0.5,0,0.5,1}} which counts the proportion of cell interfaces manually labelled as p* and detected as p in the whole population of interfaces. This table is presented in Appendix 4—figure 3E.

To obtain the results presented in Figure 3—figure supplement 1G, we chose to consider as ‘Correct’ predictions not only the cases where p*=p, but also when expert and algorithm agree on direction, that is when p* and p have the same sign, independently of confidence. If we sum up those categories, we end up with more than 72% of cell interfaces being predicted correctly (green squares in Appendix 4—figure 3E). Conversely, we considered predictions to be ‘Opposite’ to the ground truth when p* and p have opposite signs (red squares in Appendix 4—figure 3E) which only accounts for less than 10% of the interfaces. These 10% of incorrectly predicted interfaces are a problem that may come from an increased sensitivity to image noise with a coarser resolution, but their effect might be mitigated by the integration of interface polarities at cell level (see below), especially since they tend to concern interfaces of smaller area (Appendix 4—figure 3F).

The rest of the interfaces was considered as ‘Uncertain’ as they either consist of interfaces for which the experts determined a polarity when the automatic method did not lead to statistical significance (almost 13% of the interfaces) or the opposite case of detecting a polarity when the experts were consistently uncertain (the remaining 5%). In the first case, the interfaces will have a mild effect on the resulting cell level vectors as they will contribute as a null vector.

This mitigated effect of interface-level polarity errors when we integrate polarity at cell level was evaluated by computing a less quantitative version of the cell polarity vectors than the one presented in the next section. The polarity vector of the cell labelled c is noted 𝐏𝐈𝐍c and is computed as an area-weighted average of interface normals multiplied by their polarity value:

(76) 𝐏𝐈𝐍c=c𝒜c1Ac,cpcc𝐧ccc𝒜c1Ac,c

This allows to compare the directions of the resulting cell vectors. In Appendix 4—figure 3G–J, we can see how cell interface polarities that were not correctly predicted (orange circles in Appendix 4—figure 3H) lead to cell polarity vectors that have very close directions to those resulting from the ground truth polarities. This preservation of cell polarity vector was studied more extensively on the 83 cells that shared consensus interfaces (Figure 3—figure supplement 1H). Appendix 4—figure 3K shows that more than 50% of the considered cells have their polarity vector oriented in the exact same direction than the ground truth while than more than 90% were found in qualitatively the same direction (angle error less than 90 difference).

Appendix 4—figure 3
Evaluation of cell interface PIN1 polarities using manually annotated SRRF images.

(A) The super-resolution radial fluctuations method produces a 2D image slice with a greatly improved pixel resolution and a less spread-out signal. (B) Using the ImageJ line tool, polarities were manually estimated for each pair of neighbor cells by looking at the relative position of the peaks of PIN1 and PI signals along lines locally orthogonal to the cell interface. In the displayed case, the ground truth polarity value would be -1 as PIN1 lies clearly and consistently to the right of PI. (C) The interfaces on which both experts were in consensus form a set of polarity values that we use as a ground truth to evaluate the results of the automatic method. (D) Comparison of cell interface polarities assigned by the first expert taken as reference (y-axis) with the ones assigned by the second expert (x-axis). Each cell corresponds to the percentage of annotated interfaces (N = 110) for which expert one assigned the polarity y and expert two the polarity x. Agreement between experts is represented in green cells (85%), disagreement in dark red cells (4%) and uncertain cases in orange cells (11%). Empty cells correspond to a percentage 0%. (E). Similarly displayed comparison of ground truth polarities taken as reference (y-axis) and detected polarities. More than 72% of interfaces have ‘Correct’ predicted polarities (green cells) and less than 10% have ‘Opposite’ polarities (dark red cells). 18% of cells have an ‘Uncertain’ polarity decision (orange cells). (F). The interfaces where errors were made appear to be generally smaller ones. (G–J). Local errors between the ground truth interface polarities (G) and automatically detected ones (H) as those highlighted in orange circles do not necessarily translate into errors between the corresponding polarity vectors integrated at the level of the cells (I and J). (K). Cumulative histogram of cell-level angle errors, showing a very high angular consistency at cell level despite the proportion of interfaces with opposite polarities.

Computation of cell polarity vectors

The local cell interface 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 c as a 3D vector 𝐏𝐈𝐍c that takes into account all the anticlinal cell interfaces surrounding the considered cell. Each interface 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:

(77) 𝐏𝐈𝐍c=1c𝒜c1Ac,cc𝒜c1Ac,c(|PINcc-PINcc|pcc𝐧cc)

Polarity vector alignment

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. The processes of Rigid time registration and Sequence SAM reference frame determination allow to transform rigidly the images from their original frame to the common reference SAM frame *. Therefore, to transform the PIN polarity data into the same referential, we apply this transform to the cell barycenters Pc and to the cell interface centers Pc,c and apply the rotation component of the transformation to the PIN cell polarity vectors PINc,c.

Polarity vector maps

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 quantification of PIN1 signal at the level of cell interfaces and the estimation of interface areas:

(78) PIN^(𝐱)=1c1c𝒜c1η(𝐱-Pc,c*)Ac,cc1c𝒜c1η(𝐱-Pc,c*)Ac,cPINc,c

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:

(79) 𝐏𝐈𝐍^(𝐱)=1c1η(𝐱-Pc*)c1η(𝐱-Pc*)𝐏𝐈𝐍c*

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 𝐏𝐈𝐍^(x,y)=(PIN^x,PIN^y)(x,y). This operation gives us a scalar map (as those displayed in Figure 3—figure supplement 2, third column) where negative values correspond to areas of local convergence of the PIN directions:

(80) div𝐏𝐈𝐍^(x,y)=𝐏𝐈𝐍^xx(x,y)+𝐏𝐈𝐍^yy(x,y)

We use this computed map to approximate the values of div𝐏𝐈𝐍^ at punctual locations such as nuclei points Pn or primordia extremal points (ri-,θi-) or (ri+,θi+) by averaging the values at the four closest grid coordinates (Figure 3J). Indeed, it would be very complex to estimate this quantity on a discrete set of points from the 𝐏𝐈𝐍c* vectors, on which the divergence operator would be non-trivial.

Implementation details

Image processing

Histogram normalization was performed using the implementation provideed in the Scikit-Image library (van der Walt et al., 2014), and cell barycenters in the segmented images were computed using the SciPy library (Virtanen et al., 2020), made available through the timagetk (Tissue Image ToolKit) image processing library.

Meshing algorithms

We used the implementations provided in the VTK (Visualization ToolKit) library (Schroeder et al., 2006) for the Marching Cubes algorithm, the mesh smoothing and the quadric decimation (vtkImageMarchingCubes, vtkWindowedSincPolyDataFilter, vtkQuadricDecimation). A version of the isotropic remeshing algorithm has been implemented in the cellcomplex library.

PIN polarity quantification

We have developed specific tools for the quantification of PIN1 polarities, and they are made avaliable in the tissue_paredes library. We used the implementation of Wilcoxon test from the SciPy library (Virtanen et al., 2020).

PIN polarity maps

Continuous maps and vector maps were generated using the functionns implemented in the tissue_nukem_3d Python library. The divergence maps were estimated by applying a 1-D Gaussian derivative filter (σ=3.75μm) to each component of a vector map in the corresponding dimension, using functions implemented in the SciPy library (Virtanen et al., 2020).

Appendix 5

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 plastochron variability make the observed systems close to a steady developing regular phyllotactic system with a divergence angle α=137.5±6.7 and a plastochron T=12h±2h (Figure 1—figure supplement 2D–G, Appendix 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 S must follow the properties of a spatio-temporally periodic function of spatial period (0,α,0) and temporal period T:

(81) (rθz)+×[-π,π]×,t𝒯,S((rθz),t+T)=S((rθz)+(0α0),t)

This helped us consider signal dynamics on durations that largely overpass the observation range of 10 to 14 hr, by applying successive rotations of the meristems to simulate the passing of time. Notably, an aligned SAM observed at t=0h rotated of 137.5° clockwise has been shown to be the best next frame to the same aligned SAM observed at T=10h (Figure 1G, Figure 1—figure supplement 2H). More generally, any primordium of stage p visible at time t can be used to infer information at time t+pT. This means that we were able to reconstruct trajectories of signals at a given 2D projected position (r,θ) over a duration of up to nine plastochrons (∼ 100 hours) from observations spanning only one, but with nine visible primordia -3 to 5. This was achieved only by interpolating rotated aligned information (Video 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 p, then they can be positioned on the same developmental axis at a time t+pT 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 hr (Figure 2G). But to achieve the mentioned 100 hr 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 (Appendix 2) in which the azimuthal component is in average close to 0, with a limited amplitude compared to the radial component (Appendix 5—figure 1A). We note 𝐯n the speed of the nucleus labelled n in the 2D polar SAM reference frame *:

(82) n1,𝐯n=vr,n𝐫(Pn)+vθ,n𝜽(Pn)

where 𝐫(Pn) and 𝜽(Pn) are the local normal and tangential unitary vectors at Pn . In the following, we will then assume a pure radial motion of cells in the L1 of the SAM, that is that:

(83) n1,vθ,n0

We compute the local radial cellular speed as a 2D continuous map on the L1 (Appendix 5—figure 1B-C), using the same parameters as for the other cellular signals: (r,θ)2,

(84) vr^(r,θ)=1n1η((r,θ)-Pn)n1η((r,θ)-Pn)vr,n.

Defined this way, local radial cellular speed is a tissue-level information, that is not attached to a cell but to a spatial position (r,θ). Therefore it has a spatio-temporal periodicity property and we can write:

(85) vr^(r,θ,t+T)=vr^(r,θ+α,t)

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 {t0<<tn<T}. Let us consider a cell with an initial position P(t0)=(r(t0),θ0), setting r(t0)=r0. Using acquisitions at t0, it is possible to estimate vr^(r0,θ0,t0), which we use to estimate P(t1)=(r(t1),θ0) assuming a linear motion between t0 and t1:

(86) r(t1)=r0+(t1-t0)vr^(r0,θ0,t0)

More generally, with observations at ti-1, we derive P(ti) with:

(87) r(ti)=r(ti-1)+(ti-ti-1)vr^(r(ti-1),θ0,ti-1)

We perform this progression until we compute the last position P(tn) 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 plastochron, and estimate P(t0+T) with:

(88) r(t0+T)=r(tn-1)+(T-(tn-1-t0))vr^(r(tn-1),θ0,tn-1)

To proceed further, we would like to progress in time and estimate the cell position at t1+T. This is where we use the spatio-temporal periodicity property to derive that:

(89) r(t1+T)=r(t0+T)+(t1t0))vr^(r(t0+T),θ0,t0+T)=r(t0+T)+(t1t0))vr^(r(t0+T),θ0+α,t0)

We are therefore able to estimate this next position from observations at t0 simply by rotating the radial speed map (Appendix 5—figure 1D). Then iteratively it is possible to go further in time to tn+T, extrapolate motion to t0+2T, apply again spatio-temporal periodicity to reach t1+2T, and so on. Finally we obtain the two following general equations:

(90) i[[1,n]],r(ti+pT)=r(ti-1+pT)+(ti-ti-1))vr^(r(ti-1+pT),θ0+pα,ti-1)
(91) r(t0+(p+1)T)=r(tn-1+pT)+(T-(tn-1-t0))vr^(r(tn-1+pT),θ0+pα,tn-1)

These are valid for positive integer values of p until reaching the maximal extent in the considered data (i.e. when the map used to estimate local radial speed becomes ill-defined), which determines a value pmax. This defines a radial trajectory in the 2D space that reflects the local motion of cells along several plastochrons (Appendix 5—figure 1E). Note that a rigorously identical approach can be used to go backwards in time with negative values of p until t0+pminT, using motion vectors computed using inverse image deformation (Appendix 2).

In the end, from a single initial position (r0,θ0) we obtain a discrete radial trajectory:

(92) {(r(ti+pT),θ0)i[[0,n]],p[[pmin,pmax]]}

To monitor a cellular process over long time courses, the objective would be to estimate the value of the signal S along this spatio-temporal trajectory, that is:

(93) {S^(r(ti+pT),θ0,ti+pT)i[[0,n]],p[[pmin,pmax]]}

Using the spatio-temporal periodicity property of S, this translates into:

(94) {S^(r(ti+pT),θ0+pα,ti)i[[0,n]],p[[pmin,pmax]]}

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 (Appendix 5—figure 1F-H). This is the approach we use to reconstruct cell-level 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).

Appendix 5—figure 1
Using spatio-temporal periodicity to extrapolate long-term cellular motions.

(A). Cellular motion on the L1 is essentially a radial motion towards the periphery of the SAM that accelerates as cells exit the CZ. (B–C). Average 2D maps of L1 local radial cellular speed, computed from image deformations between observations at t = 0 hr and t = 5 hr (B) and t = 5 hr and t = 10 hr (C) (N=21 SAMs). (D). Spatio-temporal periodicity allows estimating cellular motion on several plastochrons by successive rotations: motion at t0 can be used to estimate position at t1, which can be extrapolated to t0+T. By periodicity, radial motion a t0+T is equal to radial motion at t0 rotated by one divergence angle α, which gives the position at t1+T, and so on. (E). The iteration of this process in time allows the reconstruction of long-term radial cellular trajectories that correspond to the average motion of cells over the population over nearly 100 hr. (F–H). These trajectories are used to define spatial domains that reflect cellular motion over time and allow the study of cell-level processes over long durations.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
    Proceedings of the 2004 eurographics/ACM SIGGRAPH symposium on geometry processing- SGP ’04
    1. M Botsch
    2. L Kobbelt
    (2004a)
    Proceedings of the 2004 Eurographic.
  9. 9
    A remeshing approach to multiresolution modeling
    1. M Botsch
    2. L Kobbelt
    (2004b)
    Proceeding of the Symposium on Geomentry Processing. pp. 185–192.
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
    Proceedings of the 24th annual conference on computer graphics and interactive techniques
    1. M Garland
    2. PS Heckbert
    (1997a)
    SIGGRAPH’ 97 ACM Press.
  28. 28
    Surface simplification using quadric error metrics 
    1. M Garland
    2. PS Heckbert
    (1997b)
    Proceedings of the 24th Annual Conference on Computer Graphics and Interactive Techniques SIGGRAPH. pp. 209–216.
  29. 29
  30. 30
    czifile Laboratory for fluorescence dynamics University of California
    1. C Gohlke
    (2012)
    czifile Laboratory for fluorescence dynamics University of California, http://www.lfd.uci.edu/.
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
    Marching cubes: a high resolution 3D surface construction algorithm
    1. WE Lorensen
    2. HE Cline
    (1987)
    Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques SIGGRAPH. pp. 163–169.
  42. 42
  43. 43
  44. 44
    Spatio-Temporal registration of 3D microscopy image sequences of Arabidopsis floral meristems
    1. G Michelin
    2. Y Refahi
    3. R Wightman
    4. H Jonsson
    5. J Traas
    6. C Godin
    7. G Malandain
    (2016)
    Proceedings - International Symposium on Biomedical Imaging. pp. 1127–1130.
  45. 45
  46. 46
  47. 47
    Medical Image Computing and Computer-Assited Intervention
    1. S Ourselin
    2. A Roche
    3. S Prima
    4. N Ayache
    (2000)
    557–566, Block matching: A General Framework to Improve Robustness of Rigid Registration of Medical Images, Medical Image Computing and Computer-Assited Intervention, 1935, Springer, 10.1007/978-3-540-40899-4_57.
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
    The Visualization Toolkit–An Object-Oriented Approach to 3D Graphics, 77
    1. W Schroeder
    2. K Martin
    3. B Lorensen
    (2006)
    Kitware, Inc.
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
  71. 71
  72. 72
  73. 73
  74. 74
  75. 75
  76. 76
    PIN-FORMED 1 regulates cell fate at the periphery of the shoot apical meristem
    1. T Vernoux
    2. J Kronenberger
    3. O Grandjean
    4. P Laufs
    5. J Traas
    (2000)
    Development 127:5157–5165.
  77. 77
  78. 78
  79. 79
  80. 80
  81. 81
  82. 82

Decision letter

  1. Jürgen Kleine-Vehn
    Reviewing Editor; University of Natural Resources and Life Sciences, Austria
  2. Christian S Hardtke
    Senior Editor; University of Lausanne, Switzerland

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

Your impressive quantitative time lapse approach provides novel insight into the dynamics of auxin distribution in shoot apical meristems. This is an important study for the field of auxin biology and certainly defines a new standard.

Decision letter after peer review:

Thank you for submitting your work entitled "Temporal integration of auxin information for the regulation of patterning" for further consideration by eLife. Your article has been evaluated by Christian Hardtke (Senior Editor) and Reviewing Editor Jürgen Kleine-Vehn.

The reviewers were highly positive about your work. The reviewers do not request additional experiments, but ask you to revise your manuscript for style and clarity. eLife is certainly a forum that allows you to highlight the novelty of your findings, but at the same time also to critically discuss your data and experimental shortcomings. We are very much looking forward to receiving a revised version of your manuscript, which integrates the constructive comments of the reviewers (see below).

Reviewer #1:

Auxin is a central determinant of pattern formation in plants. Here Galvan-Ampudia et al. use a ratiometric DII/mDII reporter to address how spatio-temporal pattern of auxin emerges. It is not very surprising that spatial auxin distribution (as visualized by DII) in the shoot apical meristem follows a precise reiterative pattern, since this correlates with (precise) organ formation. However, the authors here visualize the emergence of this pattern with a high temporal and cellular resolution. However, the manuscript is at places hard to follow and could benefit from better explanations of the methodology and thorough discussion of the related findings. A central claim of this study is that auxin maxima travel through the meristem, which cannot be explained by growth, but depends on intercellular transport. This again confirms previous assumptions, but it would be nice to better describe the underlying data. I am for example wondering how the confinement or enlargement of an auxin maxima would affect the calculations/predicted speed of auxin waves. In the given example (Figure 2A-C), the auxin maxima/DII minima seems to be in time confined to less cells. Do the authors define such a confinement as traveling?

Especially in the central zone, DII dynamics do not correlate with DR5 dynamics confirming the anticipated auxin signaling repression in this zone. The authors experimentally address its nature by exogenous application of auxin. I am a bit worried about the very strong conclusions related to these experiments. Considering that the tissue layout and cell identities differ substantially, the uptake of IAA can be highly diverse in the different regions of the SAM. Hence, I doubt that these results are easy to interpret. It is probably also challenging to obtain uptake information. Moreover, additional components (e.g. metabolism/intracellular transport) may additionally cloud the conclusions on auxin signaling mechanisms. All these aspects also affect the interpretation of the pid mutant, because it has a more regular tissue layout and unified tissue identity. The data on ETTIN is poorly discussed and I feel even misplaced. ETTIN is a non-canonical component of auxin signaling and hence may not be the best choice to discuss general aspects of auxin signaling.

Reviewer #2:

Manuscript by Galvan-Ampudia is an interesting story that demonstrates how auxin concentration is integrated with time and space to deliver new organs in the aerial section of the plant. I was particularly impressed with the combination of imaging techniques and time lapses that were used to gather quantitative data on auxin dynamics. In short, I believe it is an important study that will receive further attention. However, I would like to raise some questions that authors should address before the publication:

1) Results paragraph two:

It is not clear to me what authors mean by minimal information loss when they align SAM images (percentage of error?). This should be clarified and explained better for a general reader.

2) A dynamic range of DII-VENUS is not ideal. In fact, this reporter can only discern between relatively low auxin concentrations – we can actually see some signal. Therefore, it is difficult to distinguish between intermediate or high auxin concentration. Traveling auxin waves would be expected to carry rather medium or high auxin concentrations. Could the authors elaborate more on this matter?

3) Results paragraph four: Has to be clarified – this sentence is difficult to understand

4) I do not understand this sentence "reconstruction of several initiation events from only one?”

5) Auxin is known to travel with the speed of ~10mm/h (Kramer et al., 2011), whereas authors observe waves moving with speeds <1µM/h. Would authors come up with a plausible explanation for this severely reduced auxin waves velocities (~ 4 orders of magnitude difference).

6) Another interesting aspect is that transcriptional responses to auxin are delayed by roughly 12h. Now, is there anything known about the mechanism that could cause such an extensive delay?

7) Auxin treatments presented in the manuscript seem to involve unusually high if not extreme auxin concentration (mM range). Are those concentrations physiologically relevant? Could those effects observed be merely an artifact that is caused by overdosing of auxin?

8) Subsection “Spatio-temporal control of auxin efflux and biosynthesis” paragraph two: Interesting concept presented there. Would then PINs align towards or alternatively from an auxin source (production sites) based on new data presented in the manuscript.

9) Subsection “The role of time in transcriptional responses to auxin”, final sentence in first paragraph: very unclear sentence that needs revising.

Reviewer #3:

This manuscript by Galvan-Ampudia et al. reports an intensive study that investigated spatio-temporal patterning of auxin activity in the shoot meristem. The data reveal an intriguing dynamic over space and time that, to my knowledge, is captured here for the first time. There are various interesting observations here, for instance the observed partial temporal uncoupling between auxin levels and auxin response, which provides important information for judging auxin response in different contexts. Overall, the effort the authors made to monitor auxin dynamics quantitatively in a growing meristem are impressive. There are some limits one could discuss, for example the methodological approach to determine PIN polarity, but in any case, in my opinion the methods are state-of-the-art and thus as good as it gets for now. The manuscript is well written and documented, although I must add that I cannot judge the modeling part competently. It also shows that this paper has already undergone two rounds of revision elsewhere, it is very much honed and as far as I am concerned ready to publish.

https://doi.org/10.7554/eLife.55832.sa1

Author response

Reviewer #1:

Auxin is a central determinant of pattern formation in plants. Here Galvan-Ampudia et al. use a ratiometric DII/mDII reporter to address how spatio-temporal pattern of auxin emerges. It is not very surprising that spatial auxin distribution (as visualized by DII) in the shoot apical meristem follows a precise reiterative pattern, since this correlates with (precise) organ formation. However, the authors here visualize the emergence of this pattern with a high temporal and cellular resolution. However, the manuscript is at places hard to follow and could benefit from better explanations of the methodology and thorough discussion of the related findings.

We apologize for the lack of clarity in some parts of the manuscript. This is largely due to the original choice of a short format. We have followed the comment of the reviewer to address this (see below).

A central claim of this study is that auxin maxima travel through the meristem, which cannot be explained by growth, but depends on intercellular transport. This again confirms previous assumptions, but it would be nice to better describe the underlying data. I am for example wondering how the confinement or enlargement of an auxin maxima would affect the calculations/predicted speed of auxin waves. In the given example (Figure 2A-C), the auxin maxima/DII minima seems to be in time confined to less cells. Do the authors define such a confinement as traveling?

We have measured the motion of the centers of auxin local maxima, independently of the spatial extent of these zones. We have not analyzed how the extent of auxin maxima varies over time, as auxin distribution show dynamic and complex changes over time. In the case of Figure 2A-C, we highlighted cells (green circle) that were around the center of the accumulation zone at time 0 (and not the extent of the auxin maximum zone) and end up depleted in auxin at time 10 (white circles). It can also be seen that the cells on the distal edge of the accumulation zone at time 0 (red circles) at the same time undergo a progressive increase of their auxin level. In other words they follow opposite trends in auxin. Our data thus indicate that these changes are not simply due to a confinement of high auxin concentrations to fewer cells but rather to a spatial shift in the auxin distribution pattern relatively to the cellular canvas of the tissue. We have clarified this in the main text, modified Figure 2A-C and amended the legend of Figure 2.

But it is true also that the maximum is confined to fewer cells with the progressive establishment of the auxin minimum at the boundary of the CZ. We have thus added this in the main text.

Especially in the central zone, DII dynamics do not correlate with DR5 dynamics confirming the anticipated auxin signaling repression in this zone. The authors experimentally address its nature by exogenous application of auxin. I am a bit worried about the very strong conclusions related to these experiments. Considering that the tissue layout and cell identities differ substantially, the uptake of IAA can be highly diverse in the different regions of the SAM. Hence, I doubt that these results are easy to interpret. It is probably also challenging to obtain uptake information. Moreover, additional components (e.g. metabolism/intracellular transport) may additionally cloud the conclusions on auxin signaling mechanisms. All these aspects also affect the interpretation of the pid mutant, because it has a more regular tissue layout and unified tissue identity.

We agree with the reviewer that differential uptake could be an important issue. Although it is effectively challenging to obtain a precise information, having qDII in the treated meristems allows is to show that even after 30’ DII-VENUS is degraded throughout the meristem. We were showing these data only in the form of plots in Figure 5—figure supplement 1 and have now included images (Figure 5—figure supplement 1A-H) in addition to the quantifications (Figure 5—figure supplement 1I). These data suggest that the uptake could be homogeneous. For pid mutants, the fact that all cells react with the same dynamics to the treatments suggests also an uptake throughout the peripheral zone. And this is also consistent with previous publications that we have cited (10,12,20) where activation of DR5 throughout the PZ can be observed in different genetic backgrounds.

Taking this into consideration, we have modified the text to better highlight how our results support our conclusions but at the same time to point to the fact that differential auxin uptake could still bias our results.

Concerning additional components such as auxin metabolism and intracellular transport, this experiment does not demonstrate in any way that the delay that we observe in the activation of DR5 is due to mechanisms involving only auxin signaling effectors or chromatin (as supported in part by our results with ett and TSA). We completely agree that auxin metabolism and intracellular transport could be implicated in creating the delay and in integrating the auxin concentration. What is in between auxin and DR5 remains a black box in this experiment, however, it allowed us to identify the emergent properties of the different mechanisms acting in this black box. Taking into account the comment of the reviewer, we have now better discussed this and mention auxin metabolism and intracellular transport as possible mechanisms acting in the temporal integration of auxin concentration.

The data on ETTIN is poorly discussed and I feel even misplaced. ETTIN is a non-canonical component of auxin signaling and hence may not be the best choice to discuss general aspects of auxin signaling.

While it is true that ETTIN is a non-canonical ARF, the genetic data recently published by the group of Doris Wagner (Chung et al., 2019) demonstrate clearly that ETTIN is acting redundantly with ARF4 and MONOPTEROS in promoting organogenesis. In addition, our results showing a reduced DR5 expression in the ett background support this conclusion. Following the advice of the reviewer, we have now highlighted this in the text.

It is true that confirming the ETTIN results with other genetic backgrounds with perturbation of auxin signaling will be necessary in the future (we had tried this without success due to difficulties with genetic material). We have thus now indicated the need for such an analysis in the discussion at the end of the paragraph on the mechanisms that could act in time integration of the auxin signal (Discussion paragraph two).

Reviewer #2:

Manuscript by Galvan-Ampudia is an interesting story that demonstrates how auxin concentration is integrated with time and space to deliver new organs in the aerial section of the plant. I was particularly impressed with the combination of imaging techniques and time lapses that were used to gather quantitative data on auxin dynamics. In short, I believe it is an important study that will receive further attention. However, I would like to raise some questions that authors should address before the publication:

1) Results paragraph two:

It is not clear to me what authors mean by minimal information loss when they align SAM images (percentage of error?). This should be clarified and explained better for a general reader.

We have rephrased this in the following way:

“All images could be superimposed preserving the spatial distribution of auxin maxima and minima “

2) A dynamic range of DII-VENUS is not ideal. In fact, this reporter can only discern between relatively low auxin concentrations – we can actually see some signal. Therefore, it is difficult to distinguish between intermediate or high auxin concentration. Traveling auxin waves would be expected to carry rather medium or high auxin concentrations. Could the authors elaborate more on this matter?

We are not entirely sure that we follow the reasoning of the reviewer here. Our quantification demonstrates that we are able to detect blue and yellow signal in most cells of the PZ in the SAM. In the PZ, we only get close to saturation at P0 i.e. the absolute maximum of auxin in the meristem (absolute minimum of VENUS fluorescence). 1-DII-VENUS/TagBFP is at 1 when it saturates and one can see for example Figure 4G that most of the values are below 1. qDII has thus a dynamic range that detects the concentration present in the PZ and thus to do the analysis that we report in this manuscript. To make that sure to readers we have amended the text and state clearly that the dynamic range of qDII provide information without saturation in the vast majority of cells in the PZ.

3) Results paragraph four: Has to be clarified – this sentence is difficult to understand

We have clarified the sentence in the following way: “We then considered the temporal changes in auxin distribution by using time-lapse images over one plastochron, which corresponds to the period of this rhythmic system”

4) I do not understand this sentence "reconstruction of several initiation events from only one?”

We apologize but we could not find this phrase in the main text of the manuscript and were thus not able to provide an answer to this concern.

5) Auxin is known to travel with the speed of ~10mm/h (Kramer et al., 2011), whereas authors observe waves moving with speeds <1µM/h. Would authors come up with a plausible explanation for this severely reduced auxin waves velocities (~ 4 orders of magnitude difference).

The speed compiled and reported in Kramer et al., 2011 come from experiments done on tissue segments (root, stem) and very likely involve the vascular tissue. So they give an idea of the speed of transport on long distances going at least in part through specialized tissues, which could be part of the explanation.

More importantly, the speed that we estimate is not the speed at which auxin itself travels between cells but a motion of the distribution pattern. The spatial pattern of the auxin distribution can be seen as the current “steady state” of a system where, among other processes, cells are constantly exchanging auxin through transport at a much higher rate. The dynamic reconfiguration of this system (notably through the changes in orientations of the PIN1 network) leads to changes in its steady state, which takes place over much longer time frames. The motion of auxin accumulation zones corresponds to the spatial dynamics of this slowly varying steady state, hence the smaller speeds.

For these different reasons, we simply think that these are two type of data that cannot be compared. It seemed thus difficult to include a specific text on this point in the manuscript. We have simply stated on several occasion that velocity/movement/wave are only “apparent”.

6) Another interesting aspect is that transcriptional responses to auxin are delayed by roughly 12h. Now, is there anything known about the mechanism that could cause such an extensive delay?

It is the first time that such a delay is demonstrated in response to auxin or to our knowledge to any other developmental signals in plant. It is beyond the scope of the manuscript to elucidate the mechanism beyond this delay and at this stage we can only discuss/speculate about the mechanisms at play. We had and are still discussing a possible mechanism that could be based on the regulation of the acetylation state of auxin target genes. Regulation of chromatin state is a mechanism involved in time integration in animal and plants (Nahmad et al., 2011; we slightly amended the text to make that clearer). This is also pertinent in light of a recent publication by some of the authors (Ma et al., 2019) that we discuss and that show that WUS represses auxin signaling genes and auxin target genes in the CZ through histone deacetylation. The delay could then be due to the time taken to acetylate the auxin target loci. To address a comment of reviewer #1, we are also now suggesting that mechanisms regulating intracellular auxin concentrations could also regulate the delay. This can be found in the Discussion paragraph two.

7) Auxin treatments presented in the manuscript seem to involve unusually high if not extreme auxin concentration (mM range). Are those concentrations physiologically relevant? Could those effects observed be merely an artifact that is caused by overdosing of auxin?

We have used concentrations of exogenous auxin that have already been shown to be physiologically relevant for the SAM. This concentration would be unusually high for roots but they are not for the SAM where concentrations up 30 mM have been used to induce organogenesis or gene expression at the SAM (Reinhardt et al., 2000; Heisler et al., 2005, also Sassi et al. Curr Biol 2014). So the effects that we see and for which we use a concentration as low as 0.1 mM are very unlikely to be caused by overdosing of auxin.

We have thus simply amended the text in the following way: “To test these scenarios, we treated SAMs with auxin for different periods using physiologically relevant concentrations (11,12)”

8) Subsection “Spatio-temporal control of auxin efflux and biosynthesis” paragraph two: Interesting concept presented there. Would then PINs align towards or alternatively from an auxin source (production sites) based on new data presented in the manuscript.

Effectively our data suggest a situation where PIN1 proteins allow for pumping auxin away from the auxin production sites. We have added extra text stating this.

We would like to point at the same time that PIN1 polarities allows pumping auxin toward auxin maxim at the center and organ initia. So clearly further work is needed to understand what drives the organization of the polarities at the SAM scale.

9) Subsection “The role of time in transcriptional responses to auxin”, final sentence in first paragraph: very unclear sentence that needs revising.

We apologize. This was simply a typo: “at” changed for “is” to clarify the sentence.

Reviewer #3:

This manuscript by Galvan-Ampudia et al. reports an intensive study that investigated spatio-temporal patterning of auxin activity in the shoot meristem. The data reveal an intriguing dynamic over space and time that, to my knowledge, is captured here for the first time. There are various interesting observations here, for instance the observed partial temporal uncoupling between auxin levels and auxin response, which provides important information for judging auxin response in different contexts. Overall, the effort the authors made to monitor auxin dynamics quantitatively in a growing meristem are impressive. There are some limits one could discuss, for example the methodological approach to determine PIN polarity, but in any case, in my opinion the methods are state-of-the-art and thus as good as it gets for now. The manuscript is well written and documented, although I must add that I cannot judge the modeling part competently. It also shows that this paper has already undergone two rounds of revision elsewhere, it is very much honed and as far as I am concerned ready to publish.

We thank the reviewer for the very positive evaluation of the manuscript. We have taken into account the comment on PIN1 polarity. Concerning our approach and to make clear that there some limits to it, we now state the following:

“This quantitative evaluation (Figure 3D-G, Figure 3—figure supplement 1 and Appendix 4) validates the robustness of our method, showing that, in spite of a coarse image resolution, a vast majority of cellular polarity directions are consistent with super resolution imaging techniques. Our approach is therefore particularly suitable for monitoring global trends at the scale of a tissue”.

[Editors' note: we include below the reviews that the authors received from another journal, along with the authors’ responses.]

First round of Review

Reviewer 1:

In plants, the hormone Auxin serves as a major regulator of development, growth and physiological response. Here, looking at organ primordium formation in the Arabidopsis shoot meristem (SAM), the authors address the question of whether the helical arrangement of organs reflects positional information provided by variations of auxin levels in time in addition to (as was previously described) in space. This is an exciting and provocative story for several reasons: (1) it challenges ideas about the primary role of the PINs in auxin-mediated development; (2) it describes a new time-dependent patterning phenomenon and (3) suggests a mechanism for how “memory” of prior auxin can lead to a delayed response. The imaging data are of high quality and these and the modeling data are, in general, presented in a logical and aesthetically pleasing way.

We appreciate this very positive comment on our work.

My comments related to the 3 major points of this story:

1) This is clearly challenging PIN-based SAM models and I expect the experts in those fields/proteins can evaluate the subtle differences in behaviors better than I, and also say whether this paper fairly characterizes previous data and interpretations. I will restrict myself to evaluating the data presented here and the authors conclusions. As an outsider, whether they disprove specific PIN hypotheses isn’t all that important to me-if that’s all this paper did, I wouldn’t find that alone to be sufficiently compelling, and it they failed to disprove an old hypothesis I could see ways that this would be suitable for publication.

We thank the reviewer for his appreciation of this part of our work. We agree that the point of our manuscript is not to challenge the PIN1-based SAM models. The reason for exploring the distribution of PIN1 polarity is rather linked to the unexpected dynamics of auxin distribution we measured and that cannot be easily explained by the existing models. We would like to point out here that we have now further validated our approach for quantifying PIN1 polarities using super resolution microscopy. A comparison between the polarities computed from a 3D confocal image and the one determined using single plane super resolution image is now shown on Figure 3B-G. We have also added a new Figure 3—figure supplement 1 with further images and a quantification of the differences between polarities determined with the two approaches. The quantification shows that more than 80% of PIN1 cellular polarities match almost perfectly (deviation less than 30°) when comparing the two approaches (Figure 3—figure supplement 1 and corresponding text).

Note also that we now provide genetic data that, in addition to published work, suggest that auxin biosynthesis from YUC1 and YUC4 in flowers is important for phyllotaxis. Taken all together, our findings provide a tissue-scale vision of the organization of PIN1 polarities and auxin biosynthesis that is compatible with the auxin distribution we observe.

2) This is the strongest portion of the manuscript. It is based on the generation of new reporters (a ratiometric DII-FP) and imaging/computational innovations that showed that auxin distribution (as measured by this reporter) was much more stereotyped than one might have naively expected given that plant cell divisions in the SAM are not predetermined (i.e., the patterning is not via a C. elegans-like invariant lineage). Still, following 21 SAMs over a 10 hr interval, it was possible to see that the same reporter behavior almost to cell level resolution. This consistency allowed the authors to explore a time dimension of auxin-mediated patterning, and to test previous models of how auxin maxima are associated with cellular growth. In general, I found the experimental data reported here (Figures 1-4) convincing.

We thank the reviewer for this positive evaluation and agree that this is a crucial part of our analysis.

3) It is convincing that DR5 and DII-Venus have complex, non-linear relationships. It is a really exciting possibility that there is something happening in the SAM that enforces a delay in response. However, non-linear relationships could also be due to less direct or interesting mechanisms than the authors suggest, and of course “inputs” and “responses” are based on imperfect sensors of auxin levels and auxin response. Still, a stronger argument for auxin priming responses in cells could be obtained by modifying the auxin treatments. Of the two possibilities presented, the first one is never really tested. For the second, it’s not really clear that its time and not overall amount of auxin (amount integrated over time). I compare this to experiments with light responses where people are careful to include controls the test the total # of photons received to ensure that a particular response is due to duration of light exposure and not just total amount. Could you take advantage of the regularity of the SAM to compare short pulses of high auxin vs. longer pulses of lower amounts to test the differential effects on different cell populations?

We agree with the reviewer that analyzing the effect of auxin treatments with different concentrations and durations on the induction of transcription was an excellent way to consolidate this essential part of our manuscript. We have now included the results of new treatments with different auxin concentrations and different durations of treatments on Figure 5I. This figure now shows the effect on DR5 expression of treatments with 0.2mM, 1mM and 5 mM auxin for different durations. The effect of each treatment is shown throughout the PZ (x axis with indication of the position of primordia) as the log of the difference between DR5 expression after and before the treatment. This allows us now to show that both the time of exposure and the concentration of auxin control the activation of transcription monitored by DR5, thus supporting the conclusion that cells integrate the concentration of auxin over time and not only time.

Regarding the “first possibility” mentioned by the reviewer i.e. the possibility that the delay in activating transcription is due to intrinsic auxin-independent differences in cells, we would like to mention that we had already stated in the text that the exogenous auxin treatments on wild-type meristems addressed at least in part this question. Indeed in this scenario, the capacity to respond is linked to an intrinsic auxin-independent regulation of the competence that would create an asymmetry between different positions at the periphery of the central zone. It is very unlikely that even a 5h treatment would allow to induce transcription throughout the peripheral zone. One would expect that asymmetries in the response would still exist. Our data from Figure 5 then suggest that such asymmetries do not exist. Also, this conclusion is supported by published data (Reinhardt et al., 2003) demonstrating that in pin-shaped mutants all cells at the periphery of the meristem can respond to auxin by producing outgrowth. This suggests again that cells exiting the stem cell niche are equally competent to respond to auxin. We have now complemented these data by analyzing transcription activation by auxin (using DR5) in pinshaped inflorescence of the pinoid mutant. Organ initiation is inhibited in pinoid meristems and our results indicate that there is no pre-established auxin distribution. We also show that all cells at the periphery of the meristem can respond similarly to an exogenous auxin treatment but again activation of transcription by auxin is dependent on the duration of the treatment. We have included these new results in Figure 5J-M and Figure 5—figure supplement 2A-C.

We explain the effect of auxin treatments in wild-type and pinoid meristems in the text.

The way this section is written gives the impression that the DR5 vs. DII expression in the SAM is a unique feature of this developmental situation, Please cite the studies that show this is not the case elsewhere in the plant.

We have never claimed that a dependence on the time of exposure to auxin for DR5 activation was unique to the SAM and we are unsure why the reviewer had this feeling. Published data (such as the one published by some of the authors in Brunoud et al., 2012 and Band et al., PNAS 2012) show that DR5 might respond with no delay to endogenous or exogenous modifications of auxin concentrations. However, demonstrating that DR5 responds to auxin with no delay in the root and in other parts of the plant would require an extensive quantitative analysis that goes beyond the scope of the manuscript. We have thus highlighted this in the Discussion in the following way:

“Whether integration of auxin temporal information exists in other tissues remains to be established.”

The causal relationships between auxin, HDACs and DR5 response needs to be tested more clearly. How do you know that chromatin modifications aren’t auxin-independent events happening as cells leave the CZ? I would be more convinced by a genetic experiment that should it was an ARF-mediated effect on chromatin needed for response, not just a very broad drug treatment that potentially affects most of the chromatin.

There are few auxin signaling proteins for which a clear function has been established in the SAM. The one with the clearest implication is ARF5/MP but the mutant was shown to be insensitive to auxin in the meristem (Reinhardt et al., 2003) and thus cannot be used here. ARF3/ETTIN and ARF4 are the two other ARFs that show high expression during organ initiation (Vernoux et al., 2011). A key role for ARF3/ETTIN in organ initiation at the SAM has recently been demonstrated by a publication from the laboratory of Doris Wagner. Following the advice of this reviewer (and of another one), we have now included data demonstrating the ARF3/ETTIN gene is necessary for the integration of the auxin signal. These new data are included in Figure 5 (Figure 5N-Q) and the new Figure 5—figure supplement 2 (Figure 5—figure supplement 2D,E, and described in the text in the following way:

“The Auxin Response Factor (ARF) ETTIN (ETT/ARF3) plays an important role in promoting organogenesis in the SAM (Wu et al., 2015; Chung et al., 2019). We found that in a loss-of function ett3 mutant the expression of DR5 was restricted to only 2-3 cells at sites of organogenesis, and that a 300’ 1mM auxin treatment did not induce DR5 in the SAM (Figure 5N-Q and Figure 5—figure supplement 2D-E).”

Concerning the link to histone acetylation, we believe that it is beyond the scope of the revision of this manuscript to explore in depth whether histone acetylation levels provide the mechanism that explains time integration of auxin information. While the functional link between auxin signaling and histone acetylation is well established, understanding it further is a full research question that certainly requires years of research.

Bearing these considerations in mind, we now show the results of the pharmacological inhibition of HDACs on DR5 expression in a supplementary figure (Figure 5—figure supplement 1) and discuss them in the following after the description of the result with ARF3:

“Auxin signaling and ARF3 have been shown to act by modifying acetylation of histones (Wu et al., 2015; Chung et al., 2019 and Long et al., 2006). Pharmacological inhibition of histone deacetylases (HDACs) alone was able to trigger concomitant activation of DR5 at P0 and P-1 sites in the SAM (Figure 5—figure supplement 1I-K). These results suggest that auxin signal integration depends on a functional ARF-dependent auxin nuclear pathway.”

To further comment on this point, it is indeed likely that chromatin modifications occur both in an auxin-independent and auxin-dependent fashion. Auxin-independent changes when cells exit the stem cell niche are actually suggested by a recent work (Ma et al., 2019) showing that WUSCHEL is required to regulate the histone acetylation status at loci for a large number of auxin signaling genes and auxin targets. It is thus possible that cells exiting the stem cell niche can already experience changes in histone acetylation status due to a lower influence of WUSCHEL. We believe that, given the existing literature, this remains an important discussion point. We have thus kept the corresponding paragraph in the Discussion.

Minor notes:

Define azimuth to make more readable to non-expert audience

We added the following explanation in the text: “azimuth (angular distance)”.

Plastrocrone = plastocron?

We apologize for this. It has been corrected throughout the text.

How is this ratiometric sensor different from the other ratiometric DII sensor made by the Weijers group? Have the two been compared in the same tissue?

R2D2 from the Weijers group (Liao et al., 2015) consists in DII-VENUS and a reference fluorescent protein expressed each under their own promoter. For qDII, DII-VENUS and the TagBFP reference protein are expressed from the same promoter, to ensure an equimolar production of DII-VENUS independent of promoter activity variations.

This has been clarified in the text in the following way:

“qDII differs from previously used tools (22) in producing DII-VENUS and a non-degradable TagBFP reference stoichiometrically from a single RPS5A promoter (23, 24)(Figure S1A-H).” However, we believe that a quantitative comparison of the two type of sensors will not add much to the manuscript.

Figure 1—figure supplement 3 Croissant-shaped. Maybe crescent-shaped? Though I rather like the idea of PIN Croissants.

We apologize for this typo. During the revision of the manuscript, we have removed the corresponding term (see below also) and this has thus been corrected in this way.

Some grammar issues throughout, could be resolved through copy-editing

This has been done with the help of a native English-speaking colleague.

“Reconstruction of PIN1 polarities demonstrated that the crescent-shape often thought to indicate polarities in cells does not always correlate with polarities and can thus be sometimes misleading”. This is potentially an important technical point (how do you accurately measure cellular polarity) and it gets a bit lost here—what does “polarities in cells does not always correlate with polarities” mean? Is “tissue” missing? Given the number of people who report “PIN polarities”, I would re-emphasize the importance of accurately measuring this (and the development of such a tool to do so as reported in the supplemental material).

We agree with the reviewer that this phrase was unclear. It should have been “reconstruction of PIN1 polarities demonstrated that the crescent-shape often thought to indicate polarities in cells does not always correlate with actual polarities of the auxin flux”. In the course of the revision of the manuscript, we have removed this phrase to rather put forward the validation of the method we used to quantify PIN polarities from confocal images. As discussed earlier, we have compared the results of our method with PIN1 polarities observed using super-resolution microscopy. We agree that this was an important technical point. In the revised text, we thus highlight the fact that we have used a quantitative approach and explain how we have validated our approach.

Reviewer 2:

Unlike animals, plants reiteratively produce organs de novo from stem cells reservoirs at their growing apex, the shoot apical meristem (SAM), and much research has gone into demonstrating that auxin plays a key role in positioning the newly arising organ primordia in an orderly spatial arrangement. However, whether auxin acting as a spatial morphogen can also direct the timing of organ initiation events is unclear. Here the authors address this question by using quantitative imaging to analyze the dynamics of auxin distribution and response during organ initiation in the SAM.

This is a novel area of research that the authors tackle with invention and rigor. They for the first time accurately measure the time interval between organ initiation events at the SAM, demonstrate convincingly that auxin maxima move in waves through the SAM to regulate organ patterning, and show that auxin-induced transcription lags significantly behind auxin accumulation. In general the data are robust and sufficiently quantified.

We thank the reviewer for this very positive evaluation of our work.

However, I would also argue that in some cases the authors’ observations are not always clear and do not provide enough support for their conclusions.

It is unclear from the vaguely written Abstract what the paper is about – the authors need to be more concrete about what they show auxin does.

We have now rewritten the Abstract focusing on the main discoveries of the work i.e. the fact that auxin carries both spatial and temporal information and the demonstration that temporal integration of auxin information is necessary for activation of transcription.

The introductory text is highly technical and hard for the reader to digest. In particular, the SAM zone terminology and how the PIN1 network functions to regulate auxin distribution, promote auxin accumulation and yet also deplete auxin are briefly summarized without sufficient context.

We have modified the Introduction in order to explain better the SAM terminology. We have also added as a new Figure 1A a schematic representation of the SAM zones in order to help the reader. Concerning the role of PIN1 in regulating auxin distribution, we have better explained the role of PIN1 in controlling the orientation of the auxin fluxes and explained how the spatial organization of PIN1 polarities is thought to control auxin distribution. We believe that we now provide a sufficient context despite the constraint of the format.

In addition, because so much information is presented in such a compressed fashion in the main text, it is quite difficult for the reader to follow the data analysis and interpretation. For instance, entire complex figures, such as Figure 1—figure supplement 3 and Figure 3—figure supplement 1, are described in only a single sentence.

In the revised version of the manuscript, we have extensively edited the text (including the supplementary text and all the figure legends) to make it more accessible.

Concerning Figure 1—figure supplement 3, we effectively only summarize the main findings illustrated in this supplementary figure because of the space limitation and because our focus is primarily on the spatio-temporal dynamics of auxin maxima. However, we have now expanded the legend of Figure 1—figure supplement 3 to provide more information to the reader. We have also added longitudinal sections at primordia positions (Figure 1—figure supplement 3F-L) to illustrate better the dynamics of auxin minima. Concerning Figure 3—figure supplement 1, it has now been largely modified and simplified. To address concerns from several reviewers, we have validated the approach we used to quantify PIN polarities from confocal images. We have compared the results of our method with PIN1 polarities observed using super-resolution microscopy. A comparison between the polarities computed from a 3D confocal image and the one determined using single plane super resolution image is now shown on Figure 3B-G. We have also added a new Figure 3—figure supplement 1 with further images and a quantification of the differences between polarities determined with the two approaches. The quantification shows that more than 80% of PIN1 cellular polarities match almost perfectly (deviation less than 30°) when comparing the two approaches (Figure 3—figure supplement 1 and corresponding text). Here also, the legend provides more information about the data shown in Figure 3—figure supplement 1.

With respect to the experimentation, some of the PIN1 experiments were unclear. Figure 3A-C is described as showing PIN1 concentration over time, but only T0 data are shown.

We thank the reviewer for pointing out this error. Indeed the changes of PIN1 over time are shown in Figure 3—figure supplement 2 and not in Figure 3. We have also corrected the text accordingly notably: “Over the course of a plastochron, only limited changes of the PIN1 polarities are observed (Figure 3—figure supplement 2), …”

Also the vector field mapping experiments were poorly described and thus difficult to interpret.

We apologize for this lack of clarity. We have modified the text to describe better what we do. We do not use the term “vector fields” anymore but only talk about “vector maps” and explain in the text how we go from cell wall polarity vectors to cell polarity vectors by integrating cell wall vectors. We then explain how we generate vector maps in the following way:

“Local averaging of the cellular vectors obtained from confocal images was then used to calculate continuous PIN1 polarity vector maps in order to identify the dominant trends in auxin flux directions in the SAM (Figure 3I, and Appendix 4).”

Note that we have also reorganized the supplementary material to include an Appendix 4 that is dedicated specifically to the analysis of PIN1 polarity and provide extensive details on the methods we have used.

Figure 3 and Figure 3—figure supplement 3 do not persuasively show any YUC6 promoter activity in the SAM proper, much less in the CZ, and no evidence is provided to show that YUC- generated auxin in the primordia is important for auxin flux or organogenesis in the SAM.

We agree with the reviewer that YUC6 promoter activity was difficult to see on both figures. In the revised manuscript, we now show only YUC4 promoter activity on Figure 3 as it has the strongest activity and to limit the size of the figure. We have also modified Figure 3—figure supplement 3 to ease the visualization of the expression of YUC6. To do so, we present for YUC1,4 and 6 both the GFP signal alone (Figure 3—figure supplement 3A-C) and an overlay of the GFP signal and of the PI signal (Figure 3—figure supplement 3D-F). Concerning YUC6, a weak signal can be seen at the center of the meristem.

Moreover, we have added to the manuscript the results of an analysis of yuc1yuc4 double mutants, which showed organ positioning and morphological defects. We present these results on Figure 3L and Figure 3—figure supplement 3O-U. These genetic data, in addition to published work, support the idea that auxin biosynthesis from YUC1 and YUC4 in flowers is important for phyllotaxis. This is explained in the text in the following way: “In addition, yuc1yuc4 loss-offunction mutants show severe defects in SAM organ positioning and size (37, 38) (Figure 3L and Figure 3—figure supplement 3O-U). Taken with the organization of PIN1 polarities, these results suggest that P3-P5 are auxin production centers for the SAM that regulate phyllotaxis.”

I am also unconvinced that the data in Figure 5 demonstrate that auxin addition to the SAM for 30-120 minutes leads to a measureable increase in DR5 expression in the P0, P1 and P2 primordia because of the noise in the experimental output.

Here we made a serious mistake in the text and apologize for this. It should have been P-2, P-1 and P0 and we thank the reviewer for spotting this. In any case, we had chosen on Figure 5 a representation that was as close as possible to the original data but the noise effectively made difficult the visualization of the response for the 30 and 120’ treatments. As the noise in the logtransformed data has a symmetric distribution, we could use a statistical model to calculate the mean of the difference between DR5 expression after and before the treatment, and the 95% confidence interval of the response curves (we explain this in the Material and Method section, Quantification and statistical analysis paragraph). These data are shown on Figure 5 I. An absence of overlap of the confidence interval at a given angular position then indicates that the differences are statistically significant between treatments/control. In the case of the 30’ and 120’ treatments with 1mM auxin, significant induction can be detected at the P-1 and P-2 sites and even at the predicted angular position of P-3. A small induction is visible at P2 but there are no significant changes at P0 and P1. We explain this in the text in the following way: “In the shorter auxin treatments (30’ and 120’), the auxin transcriptional response was mainly enhanced at P-1 and P-2 and to a lesser extent at the position of the predicted P-3 i.e. where cells are already being exposed to auxin (Figure 5I).”

Note also that we have now analyzed the effect of auxin treatments on the induction of transcription using different concentrations and durations of treatment. The results are included in Figure 5I that shows the effect on DR5 expression of treatments with 0.2mM, 1mM and 5 mM auxin for different durations. This allows us now to demonstrate that both the time of exposure and the concentration of auxin control the activation of transcription monitored by DR5, thus supporting the conclusion that cells integrate the concentration of auxin over time and not simply time. In addition, we have analyzed transcription activation by auxin (using DR5) in pin-shaped inflorescence of the pinoid mutant. Organ initiation is inhibited in pinoid meristems and our results indicate that there is no pre-established auxin distribution. We also show that all cells at the periphery of the meristem can respond similarly to an exogenous auxin treatment, which provides further argument that there are no intrinsic differences in the capacity of the cells to respond to auxin. We have included these new results in Figure 5J-M and Figure 5—figure supplement 2A-C. We explain the effect of auxin treatments in wild-type and pinoid meristems in the text.

Similarly, the effect of histone deacetylase inhibition is mild at best and without additional experimental manipulation (such as the use of HDAC mutants) it seems unjustified to draw the broad conclusion that the timing of auxin signaling at a primordium depends on the chromatin acetylation status.

We agree with this reviewer and others that our analysis of the implication of histone acetylation was too preliminary to support this conclusion. We also believe that properly understanding how histone acetylation is involved in time integration would take us largely beyond what is reasonable to include in this manuscript. In the revised manuscript, we only include in the discussion the possibility that histone acetylation could provide a mechanism for time integration based on the literature as we believe that this is an important discussion point.

In link to this comment, note however that we have kept the results of the TSA treatments but rather as a way to support an implication of auxin signaling in time integration of the auxin signal. Indeed other reviewers suggested that genetics could be used to show an involvement of the auxin signaling pathway in time integration of the auxin signal. There are few auxin signaling proteins for which a clear function has been established in the SAM. The one with the clearest implication is ARF5/MP but the mutant was shown to be insensitive to auxin in the meristem (Reinhardt et al., 2003) and thus cannot be used here. ARF3/ETTIN and ARF4 are the two other ARFs that show high expression during organ initiation (Vernoux et al., 2011). A key role for ARF3/ETTIN in organ initiation at the SAM has recently been demonstrated by a publication from the laboratory of Doris Wagner (Chung et al., 2019). We have included in the revision data demonstrating the ARF3/ETTIN gene is necessary for the integration of the auxin signal. These new data are included in Figure 5 (Figure 5N-Q) and a new Figure 5—figure supplement 2 (Figure 5—figure supplement 2D,E) and described in the text in the following way:

“The Auxin Response Factor (ARF) ETTIN (ETT/ARF3) plays an important role in promoting organogenesis in the SAM (Wu et al., 2015; Chung et al., 2019). We found that in a loss-of function ett3 mutant the expression of DR5 was restricted to only 2-3 cells at sites of organogenesis, and that a 300’ 1mM auxin treatment did not induce DR5 in the SAM (Figure 5N-Q and Figure 5—figure supplement 2D-E).”

As the functional link between auxin signaling and histone acetylation is well established, we now show the results of TSA treatments on DR5 expression in a supplementary figure (Figure 5—figure supplement 1) and discuss them in the following way after the description of the result with ARF3:

“Auxin signaling and ARF3 have been shown to act by modifying acetylation of histones (Wu et al., 2015; Chung et al., 2019 and Long et al., 2006). Pharmacological inhibition of histone deacetylases (HDACs) alone was able to trigger concomitant activation of DR5 at P0 and P-1 sites in the SAM (Figure 5—figure supplement 1I-K). These results suggest that auxin signal integration depends on a functional ARF-dependent auxin nuclear pathway.”

Minor comment:

The manuscript would benefit from English language polishing, particularly to clarify word choices such as “exposition of cells to auxin” which I take to mean “exposure of cells” (Abstract and Discussion) and “diverge fluxes away from areas” where I suspect “divert fluxes away” (Discussion) is meant.

We apologize for this. We effectively meant “exposure” and “divert”. We have polished the English language as requested with the help of a native English-speaking colleague.

Reviewer 3:

This manuscript uses quantitative imaging to investigate the role of auxin in the rhythmic initiation of organs on the flanks of the shoot apical meristem. Using a ratiometric auxin sensor, the authors establish a precise spatial and temporal map of auxin accumulation in the SAM. They show that this stereotypical map is highly conserved from plant to plant, and rotates by 137 degrees (the divergence angle between consecutive organs) every 12h (the timing between the initiation of 2 consecutive organs). The authors also analyze the polarity of PIN1 polarity in the SAM and suggest a pattern of auxin fluxes different from what was previously proposed. The manuscript also describes the expression of auxin biosynthesis genes YUCCA1-11, which suggest that young primordia act as sources of auxin for the rest of the SAM, rather than auxin sinks as previously thought. Finally, the authors show that auxin response, as indicated by a DR5 reporter, is delayed compared to auxin accumulation, due to the need for cells to be exposed to auxin for a prolonged time before auxin response can take place; they suggest and propose this delay is linked to chromatin modification.

This study represents a tremendous amount of work and provides invaluable quantitative information on the role of auxin in patterning the SAM. I am particularly appreciative of the fact that the authors not only quantify auxin accumulation and response in a single representative shoot apical meristem, but do so in a large number of plants to show that the mechanisms they describe is highly conserved from plant to plant.

An important part of the data presented here is very convincing and conclusions are well supported. I believe this paper will be an important advance in our understanding of plant development.

We thank the reviewer for this very positive evaluation of our manuscript.

The link between histone deacetylation and auxin response in the SAM deserves further attention, but again, the amount of data presented in this study is impressive, and sufficient for publication.

We agree with this reviewer and others that our analysis of the implication of histone acetylation was too preliminary. We also believe that properly demonstrating how histone acetylation is involved in time integration would take us beyond what is reasonable to include in this manuscript. In this revised version, we only kept in the discussion the possibility that histone acetylation is involved in time integration based on the literature, as we believe that this is an important discussion point.

In link to this comment, note however that we have kept the results of the TSA treatments but rather as a way to support an implication of auxin signaling in time integration of the auxin signal. Indeed other reviewers suggested that genetics could be used to show an involvement of the auxin signaling pathway in time integration of the auxin signal. A key role for ARF3/ETTIN in organ initiation at the SAM has recently been demonstrated by a publication from the laboratory of Doris Wagner (Chung et al., 2019). We have included in the revision data demonstrating the ARF3/ETTIN gene is necessary for the integration of the auxin signal. These new data are included in Figure 5 (Figure 5N-Q) and a new Figure 5—figure supplement 2 (Figure 5—figure supplement 2D,E) and described in the text in the following way:

“The Auxin Response Factor (ARF) ETTIN (ETT/ARF3) plays an important role in promoting organogenesis in the SAM (Wu et al., 2015; Chung et al., 2019). We found that in a loss-of function ett3 mutant the expression of DR5 was restricted to only 2-3 cells at sites of organogenesis, and that a 300’ 1mM auxin treatment did not induce DR5 in the SAM (Figure 5N-Q and Figure 5—figure supplement 2D-E).”

As the functional link between auxin signaling and histone acetylation is well established, we now show the results TSA treatments on DR5 expression in a supplementary figure (Figure 5—figure supplement 1) and discuss them in the following after the description of the result with ARF3:

“Auxin signaling and ARF3 have been shown to act by modifying acetylation of histones (Wu et al., 2015; Chung et al., 2019 and Long et al. 2006). Pharmacological inhibition of histone deacetylases (HDACs) alone was able to trigger concomitant activation of DR5 at P0 and P-1 sites in the SAM (Figure 5—figure supplement 1I-K). These results suggest that auxin signal integration depends on a functional ARF-dependent auxin nuclear pathway.”

My main reservation is about PIN1 polarity. I do not think that confocal microscopy has a sufficient resolution to determine which side of the cell wall the PIN1 signal is located. In optimal conditions, the maximum resolution attainable with a confocal is 200nm, and the cell wall in the shoot apical meristem is barely thicker than that. I find the images and quantifications in Figure S5 B-O unconvincing at that level of resolution (also why represent the levels of PI signal as curves yet the PIN1 signal as an average value over more than a half micron?). To be fair to the authors, my reservations here are not specific to this manuscript: there are other articles (e.g. Heisler et al., 2005) presenting vector fields of PIN1 orientation/auxin flux based on images on which one cannot discern which side of the cell wall PIN1 really is, so it may very well be that earlier publications were wrong and this study is right. If anything, this study looks at PIN1 with more detail than previous studies. However, I think it still fails to solve the question of PIN1 polarity. To the best of my knowledge, the only published images clearly showing PIN1 polarity used immunolocalization (e.g. Reinhardt et al., 2003), but these are performed on sections and only provide 2D information… I think the way to solve this would be to use superresolution. Single localization microscopy (max resolution 20nm) might be hard to implement on these samples, but structured Illumination microscopy (max resolution 100nm) might allow the authors to see more clearly which side of the cell wall PIN1 is located in the L1.

We agree with the reviewer that the maximum resolution of confocal microscopes is insufficient to visibly separate the PI signal from the GFP signal on images. However, our approach uses the properties of the point spread function of light from each channel in 3D, to detect spatial shifts between the two signals throughout cell walls that would not be visible on images. This is thus a computational approach to increase the resolution and to detect polarities. We have followed the suggestion of this reviewer (and others) and used super-resolution radial fluctuations (SRRF) microscopy to validate our computational approach.

A comparison between the polarities computed from a 3D confocal image and the one determined using single plane super resolution image is now shown on Figure 3B-G. We have also added a new Figure 3—figure supplement 1 with further images and a quantification of the differences between polarities determined with the two approaches. The quantification shows that more than 80% of PIN1 cellular polarities match almost perfectly (deviation less than 30°) when comparing the two approaches (Figure 3—figure supplement 1 and corresponding text).

Note that we have also reorganized the supplementary material to include Appendix 4 that is dedicated specifically to the analysis of PIN1 polarity and provide extensive details on the methods we have used.

I am also a bit confused at the differences between PIN1 vector norm (Figure 3F), divergence map (Figure 3D) and vector shield (Figure 3B). Figure 3B shows some divergence at the center of the SAM, while Figure 3D shows lower convergence than in P-2, P-1 and P0, but Figure 3F indicates a maximum convergence at the center of the SAM. Also, Figures 3B and 3F show pattern of vector field/vector norm at P-1 pointing from the periphery to the center of the SAM, patterns that seem to correspond to divergence values between -500 and 500 in Figure 3E, yet the divergence map in Figure 3D indicates lower divergence values for P-2, P-1 and P0 than in the center of the SAM. I might not fully grasp what these panels exactly represent (the equations in supplemental data are admittedly beyond me), but this doesn’t seem fully consistent to me.

As explained above, adding the comparison between 3D quantification from confocal images and super-resolution led us to significantly modify Figure 3. We do not show anymore the data from former Figure 3B to focus our analysis on the tissue-scale organization of the polarities. Also local divergence maps (as shown in former Figure 3D) are now only shown in Figure S5.

We have also clarified how we analyze PIN1 polarities organization in the text. We do not use the term “vector fields” anymore but only talk about “vector maps” and explain in the text how we go from cell wall polarity vectors to cell polarity vectors by integrating cell wall vectors. We then explain how we generate vector maps in the following way: “Local averaging of the cellular vectors obtained from confocal images was then used to calculate continuous PIN1 polarity vector maps in order to identify the dominant trends in auxin flux directions in the SAM (Figure 3I, and Appendix 4).”

Concerning the comment of the reviewer, the PIN1 vector norm is the local length of vectors in the PIN1 vector map, which may be interpreted as the local intensity of the flux as it takes into account both PIN1 polarities and the level of the protein. The divergence map is then computed from the PIN1 vector map and thus takes into account both the directionality and the intensity information: the PIN1 polarity divergence index measures not only the organization of the PIN1 vectors but also their intensity. The higher divergence (lower convergence) in the CZ despite the visual impression on PIN1 vector maps is due to the lower PIN1 vector norm in the CZ resulting from the low expression of PIN1. We agreed that the examples displayed next to the axis of the former Figure 3E were confusing (by suggesting that divergence results only from directionality) and therefore removed them from the new Figure 3J.

Below are some more minor comments:

“Fingers-like protrusions (visible at P-2, P-1 and P0)”. P-2 is not labeled in any of the figures; please label it.

Thank you for spotting this. We have added the label to the Figures whenever we refer to this position.

Labelling the different primordia in Video S1 (and their changing stages) would also be useful.

As the stages of the primordia are changing all the time on the video, it would impossible to label them accurately. We have thus kept the video unchanged.

Line 165-166/Figure 3—figure supplement 3F: the expression of YUC6 in the CZ doesn’t show in the figure. Please show an image with only the YUC6-GFP channel and no PI.

We have modified Figure 3—figure supplement 3. We now present for YUC1, 4 and 6 both the GFP signal alone (Figure 3—figure supplement 3A-C) and an overlay of the GFP signal and of the PI signal (Figure 3—figure supplement 3D-F). Concerning YUC6, a weak signal can be seen at the center of the meristem as stated in the text.

It is not always clear whether panels in the different figures show pRPS5a::DII-Venus expression, or ratiometric values for qDII. This should be clarified.

We apologize for this lack of precision. We have now been careful to indicate in figure legends when DII-Venus, qDII ratio or Auxin (1-qDII) is presented.

Figure 1—figure supplement 2: I have trouble understanding Figure 1—figure supplement 2H. Legend indicates more/different colors than shown in the panel. Also, what is the orange dot?

We thank the reviewer for highlighting this mistake. We have corrected the legend of the figure but we only use one color now in the panel (see below).

In the previous version, Figure 1—figure supplement 2H showed blue curves for individual meristems that represented the error between the qDII map for this meristem at 10h and the map at 0h that is rotated from 180° to +180°. The red/orange curves corresponded to the same analysis but when the error was measured between qDII maps at 10h and the same 10h map that was rotated, in order to provide a reference. This analysis allows demonstrating that the error is minimized when the map at 0h was rotated from an angle close to 137°, and thus that a 137° rotation is a suitable way of approximating the temporal evolution of the system after one plastochrone. The blue dots indicated for each meristem (blue curves) the position of the minimal error. The orange dot was a mislabeling and we apologize for this.

With this comment of the reviewer, we realized that this figure panel was too complex. In the revised Figure 1—figure supplement 2H, we have thus only kept the blue curves i.e. the error between the qDII map for this meristem at 10h and the map at 0h that is rotated from -180° to +180°. The blue dots point to the error minimum when the rotation of the qDII map at 0h is close to 137°. We also provide in the figure legend a clearer explanation of what is shown in this figure panel.

Figure 1—figure supplement 3:

– The legend is insufficient and even though one can figure out the information that is not specified in the legend, it would be much easier with a more detailed legend. Please specify that the circled areas in panel A are auxin minima. More precision about what exactly the boxes in panel A would also be appreciated. I assume it represents the extent of auxin maxima?

We have expanded the legend of this figure to explain better the panel A. Concerning the boxes, they are positioned at auxin minima and maxima. They are double boxplots that represent the radial and angular variability of the values. This is also stated in the figure legend.

– I don’t understand the difference between panels B and C – they seem to show the same data with slightly different values.

We thank the reviewer for spotting this mistake. It was indeed the same data. We have corrected the figure showing the right plot in C that shows the changes in the radial position of maxima and minima.

– I find it confusing that auxin levels between auxin maxima and minima in the last P5 (last column on the right of panel 2B) seem identical: it looks like the filled box and outline box are nearly completely superimposed. It is not clear whether this is the case in panel 2C, I don’t see an outline box.

We are following the values of the maximum and minimum in P5 over the time course in Figure 1—figure supplement 3B. As stated above the corrected Figure 1—figure supplement 3C shows their radial position. The overlap for P5 is due the fact that at this stage of flower development the value of the minimum is increasing and the value of the maximum is decreasing (as can be seen in the figure), leading to a more homogeneous radial auxin distribution in the flower.

Reviewer 4:

The manuscript by Galvan-Ampudia provides convincing evidence that, in the epidermis of the inflorescence apex of Arabidopsis, cells with peak expression of a proxy of auxin concentration do not move apart because of intercalating cell growth or proliferation; rather, maxima themselves move across cells over time. The authors go on suggesting that it's the exposure of cells to auxin for a certain amount of time that defines how those cells will respond to auxin. I think this is the key and most interesting conclusion of this study

We thank the reviewer for this positive evaluation of our work and agree that these finding constitutes the key message in our manuscript.

– One, however, which I believe is not as well supported by evidence: how can the authors be sure that it's the exposure to auxin for a certain amount of time and not the total amount of auxin to which cells are exposed to over time that makes the difference? If time only were important, exposure to lower amounts of auxin over longer periods would have more dramatic effects than higher amounts over shorter periods; unfortunately, such evidence is missing from the manuscript.

We agree that consolidating this key part of our manuscript was necessary. As suggested by this reviewer as well as another one, we analyzed the effect of auxin treatments with different concentrations and durations on the induction of transcription. The new Figure 5I now shows the effect on DR5 expression of treatments with 0.2mM, 1mM and 5 mM auxin for different durations. The effect of each treatment is shown throughout the PZ (x axis with indication of the position of primordia) as the log of the difference between DR5 expression after and before the treatment. This allows us now to show that both the time of exposure and the concentration of auxin control the activation of transcription monitored by DR5, thus supporting the conclusion that cells integrate the concentration of auxin over time and not simply time. In addition, we have analyzed transcription activation by auxin (using DR5) in pin-shaped inflorescence of the pinoid mutant. Organ initiation is inhibited in pinoid meristems and our results indicate that there is no pre-established auxin distribution. We also show that all cells in the PZ of pinoid meristems can respond similarly to an exogenous auxin treatment, which provides further argument that there are no intrinsic differences in the capacity of the cells to respond to auxin. We have included these new results in Figure 5J-M and Figure 5—figure supplement 2A-C. We explain the effect of auxin treatments in wildtype and pinoid meristems in the text.

Most important, I am left wondering whether all that is biologically relevant or simply an interesting circumstance. For example, what are the effects on the development of primordia of exposure to lower amounts of auxin over longer periods and what are those of exposure to higher amounts over shorter periods? What is the biologically relevant output of those treatments?

To address this concern, we have performed daily auxin treatments of specific durations and concentrations on inflorescences from plant grown on soil in order to test how these treatments affect flower distribution on the stem and phyllotaxis. The results of these experiments support the idea that time integration of the auxin signal is biologically relevant and that its perturbation affects phyllotaxis. The corresponding data are presented on Figure 5—figure supplement 1B-D and explained in the text in the following way: “Supporting this idea, we also found that daily exogenous auxin treatments at the SAM affected phyllotaxis and that the efficiency of the treatment increased with both auxin concentration and treatment length. This was particularly evident for 30’ and 120’ treatments (Figure 5—figure supplement 1B-D). 300’ treatments were less efficient at higher auxin concentrations, possibly due to compensation mechanisms. These results suggest that temporal integration of auxin information at the SAM is essential for phyllotaxis.”

Our conclusion here is also supported by genetic data that we have added upon suggestion of the reviewer (see next point).

Moreover, one of the advantages of Arabidopsis is the outstanding genetic resources, including mutants in all the steps between perception of auxin concentrations (reported by qDII) and transcriptional responses (reported by DR5): are any of those mutants defective in their response to exposure to different amounts of auxin over different periods? And what are the biological defects associated with such defective responses?

There are few auxin signaling proteins for which a clear function has been established in the SAM. The one with the clearest implication is ARF5/MP but the mutant was shown to be insensitive to auxin in the meristem (Reinhardt et al., 2003) and thus cannot be used here. ARF3/ETTIN and ARF4 are the two other ARFs that show high expression during organ initiation (Vernoux et al., 2011). A key role for ARF3/ETTIN (together with ARF4) in organ initiation at the SAM has recently been demonstrated by a publication from the laboratory of Doris Wagner (Chung et al., 2019). In the revised manuscript, we have included data demonstrating the ARF3/ETTIN gene is necessary for the integration of the auxin signal. These new data are included in Figure 5 (Figure 5N-Q) and a new Figure 5—figure supplement 2 (Figure 5—figure supplement 2D,E) and described in the in the following way:

“The Auxin Response Factor (ARF) ETTIN (ETT/ARF3) plays an important role in promoting organogenesis in the SAM (Wu et al., 2015; Chung et al., 2019). We found that in a loss-of function ett3 mutant the expression of DR5 was restricted to only 2-3 cells at sites of organogenesis, and that a 300’ 1mM auxin treatment did not induce DR5 in the SAM (Figure 5N-Q and Figure 5—figure supplement 2D-E).”

Moreover we also added data that, together with published work (Simonini et al., 2017), show that ett-22/arf3 mutants have strong defects on phyllotaxis. These results have been included on Figure 5—figure supplement 1E-H and are explained in the text together with the results from the daily auxin treatments of specific durations and concentrations on inflorescence (see previous point). Together, these two sets of results support the conclusion that time integration of the auxin signal during organ initiation is important for phyllotaxis.

Related to this part of the manuscript, the difference between Figure 5K and J is small but visible; however, I no longer see it in the quantification (Figure 5L), which is the real data.

The difference was clearly visible and statistically significant in the quantification we presented Figure 5—figure supplement 1C in the previous version of the manuscript (now Figure 5—figure supplement 1K; see below). On Figure 5K in the previous version of the manuscript, there was effectively a script problem and a number of points were removed during the generation of the panel.

As other reviewers suggested, the link to histone acetylation was too preliminary. As properly demonstrating how histone acetylation is involved in time integration would take us beyond what is reasonable to include in this manuscript, we are not presenting the TSA treatments data in a main figure anymore. We have however kept the results of the TSA treatments but rather as a way to support the implication of auxin signaling in time integration of the auxin signal. As the functional link between auxin signaling and histone acetylation is well established, we now show the results of the TSA treatments on DR5 expression in a supplementary figure (Figure 5—figure supplement 1) and discuss them in the following way after the description of the result with ARF3:

“Auxin signaling and ARF3 have been shown to act by modifying acetylation of histones (Wu et al., 2015; Chung et al., 2019 and Long et al., 2006). Pharmacological inhibition of histone deacetylases (HDACs) alone was able to trigger concomitant activation of DR5 at P0 and P-1 sites in the SAM (Figure 5—figure supplement 1I-K). These results suggest that auxin signal integration depends on a functional ARF-dependent auxin nuclear pathway.” Note also that, in the revised manuscript, we have kept only in the discussion the possibility, based on the literature, that histone acetylation could provide a mechanism for time integration.

Still related to this part of the manuscript is the blatant absence of reference to the R2D2 sensor reported by Liao et al., 2015, to which — for all I can see — qDII is identical, except for the fluorescence protein used. This should of course be corrected.

R2D2 from Liao et al., 2015 consists in DII-VENUS and a reference fluorescent protein expressed each under their own promoter. For qDII, DII-VENUS and the TagBFP reference protein are expressed from the same promoter, to ensure an equimolar production of DII-VENUS independent of promoter activity variations.

This has been clarified in the text in the following way:

“qDII differs from previously used tools (Liao et al., 2015) in producing DII-VENUS and a non-degradable TagBFP reference stoichiometrically from a single RPS5A promoter (Wend et al., 2013; Goedhart et al., 2011)(Figure S1A-H).”

I also find interesting the part of the manuscript in which the authors suggest that claims based on qualitative visual inspection of PIN1 asymmetric distribution at the plasma membrane in the epidermis of the inflorescence apex may not always be accurate. I agree with the authors, but I don’t think their method is convincingly better: the lateral resolution of conventional light microscopy, including confocal microscopy, does not allow resolving plasma membranes of adjacent cells separated by such a thin cell wall. It would instead be necessary to compare the results of the authors' method with those of super-resolution or electron microscopy to understand the advantages and limitations of the authors' method. As mentioned, this is an interesting and meritorious part of the manuscript, but also one that is foreign to the main point of the manuscript — whether auxin concentration is only the result of transport or also that of production, for example, is not the main message of this study. And frankly I find naive that one would think that production, degradation, conjugation, etc. are irrelevant to auxin concentrations, given the strong genetic evidence that in fact they are.

We agree that a vision where only transport control the distribution of auxin is naive but the strong emphasis on the role of PIN1 in the SAM in the literature has somehow obscured the contributions of the metabolic regulation. We have attempted here to fill this gap because the spatio-temporal distribution of auxin cannot be easily explained by the dominant hypothesis present in the literature (proposing that points of convergence of PIN1 pumps lead to local accumulation of auxin). Our goal here is to provide a quantitative framework linking auxin transport and biosynthesis to auxin distribution. While we agree that it is not the main point of the manuscript, we still believe that it adds an interesting dimension to our work by providing a unique high-resolution tissue scale spatio-temporal reference dataset.

We agree with the reviewer that the resolution of confocal microscopes is insufficient to visibly separate the PI signal from the GFP signal on images. However, our approach uses the properties of the point spread function of light from each channel in 3D, to detect spatial shifts between the two signals throughout cell walls that would not be visible on images. This is thus a computational approach to increase the resolution and to detect polarities. We have followed the suggestion of this reviewer (and others) and used super-resolution radial fluctuations (SRRF) microscopy to validate our computational approach.

A comparison between the polarities computed from a 3D confocal image and the one determined using single plane super resolution image is now shown on Figure 3B-G. We have also added a new Figure 3—figure supplement 1 with further images and a quantification of the differences between polarities determined with the two approaches. The quantification shows that more than 80% of PIN1 cellular polarities match almost perfectly (deviation less than 30°) when comparing the two approaches (Figure 3—figure supplement 1 and corresponding text).

Note that we have also reorganized the supplementary material to include an Appendix 4 that is dedicated specifically to the analysis of PIN1 polarity and provide extensive details on the methods we have used.

A few minor points related to this part:

– The authors seem to use interchangeably PIN and PIN1: PIN1 is one of the 8 genes in the PIN family, so there are no PIN1 genes/proteins as the authors instead write:

It should have been PIN1 throughout the manuscript and we have corrected this.

A network of PIN-FORMED1 (PIN1) efflux carriers

– And what do the authors mean when they write that those PIN1 [sic] efflux carriers form a network? A genetic interaction network? A protein interaction network?

The organization of the PIN1 fluxes can be represented as an oriented graph whose vertices are the cell centroids and whose oriented edges between two vertices reflects the cell adjacency and the direction of the PIN1 polarity. This is the reason why we used the term network here but to avoid confusion we have removed the reference to network when talking about PIN1.

Reviewer 5:

The self-organizing properties of plant shoot meristems have fascinated researchers for decades. In the past twenty years genetic and substance-application studies have demonstrated the role plant hormones and other factors. These findings have led to complementing experimental and theoretical approaches, which link measurable parameters, like the concentrations and transport routes of the plant hormone auxin, to algorithms explaining self-organization. It is amazing that through all the thirteen years since the publication of the first of such integrated studies, many fundamental difficulties in measuring critical parameters (e.g. auxin distribution) have persisted.

It is within this context, that the submitted study fills a void as it demonstrates the successful use of new tools for a better quantification of critical parameters. The authors are introducing new markers, use them in intelligent combinations with each other and with previously characterized markers and make some stunning observations on the reproducibility of auxin distribution patterns. In particular, the use of a ratiometric auxin sensor is a huge achievement, when combined with live imaging able to show fluctuating auxin waves in the meristem at high resolution in both space and time. Moreover, by relating those patterns to auxin-response patterns at the transcriptional level, the authors raise the possibility that a cellular memory of auxin exposure has a critical role in directing a cells’ behavior.

We thank the reviewer for this very positive evaluation of our work and agree that these finding constitutes the main message in our manuscript. Note that we have consolidated our demonstration of the role of time integration of the auxin signal in organogenesis. We provide details in our response to the more specific concerns of the reviewer below.

Beyond this, the manuscript also re-assesses and re-evaluates directions of PIN1-based auxin transport routes and identifies expression sites for YUCCA genes at critical stages.

Unfortunately, these two last aspects do not add much to the overall impact of the study. The data on PIN1 polarity lead to a new picture, but are based on primary acquisition technology which is as problematic as in previous studies.

The YUCCA expression sites in developing organs, though presented as concept-changing revelations in the Abstract, may be reproducible, but could be gratuitous as there are no genetic tests in the entire paper.

Those were fair criticisms. However, we have now provided in the revised manuscript (in response of the concerns of this reviewer and others) new data showing that 1- our acquisition technology performs well and 2- that YUC expression sites are unlikely gratuitous. We provide details in our response to the more specific concerns of the reviewer below. Note also that we are providing genetic data in this revised version of the manuscript (see below).

For the four reasons listed below, I believe that the claim of the study (though surprisingly vaguely worded in title and Abstract) is too high. The authors have made great progress in measuring the dynamics of auxin distributions, and detected delay phases in transcriptional auxin responses. With a bit more cell biological work on the “hysteresis” effect, these two aspects alone would constitute a tremendous contribution.

We agree that our Abstract was too vague. We have now rewritten the Abstract focusing on the main discoveries of the work i.e., as pointed out by the reviewer, the fact that auxin carries both spatial and temporal information and the demonstration that temporal integration of auxin information is necessary for activation of transcription.

By contrast, the other, more premature parts just damage the long-term reputation of the study as a whole.

Following the reviewer advice, the part of the results that aimed at showing the role of histone acetylation in time integration of the auxin signal has been modified (see concern 3 below) and we do not claim anymore that we provide evidence for this in the manuscript.

Concerning the part on the auxin transport, while we agree that it is not the main point of the manuscript, we still believe that it adds an interesting dimension to our work by providing a high resolution tissue scale spatio-temporal reference dataset with information on auxin, distribution, its transport and its biosynthesis. It is also justified by the fact that we describe a dynamics of auxin distribution that cannot be easily explained by our current vision.

For these reasons, we have kept the analysis of auxin transport and biosynthesis in this extensively revised version of the manuscript, while at the same time focusing the logic of manuscript more clearly on the central part of the manuscript: auxin distribution dynamic characterization and the temporal integration of the auxin signal during organ initiation. As explained in more details below, all the parts of the manuscript have been consolidated.

1) Unlike for auxin concentration, the authors have no fundamentally new solution for overcoming the uncertainties in determining location, polarity and intensity of auxin transport routes and for localizing auxin sources. They end this paragraph with the revealing sentence “PIN1 in auxin distribution could also have been overestimated” by referring to not further specified “other carriers”.

We now demonstrate that we do describe a solution for a quantitative analysis of PIN1 polarity distribution. Our approach uses confocal images and the properties of the point spread function of light from each channel in 3D, to detect spatial shifts between the two signals throughout cell walls that would not be visible on images. This is thus a computational approach to increase the resolution and to detect polarities. Following suggestions from several reviewers, we used super resolution radial fluctuations (SRRF) microscopy to validate our computational approach. A comparison between the polarities computed from a 3D confocal image and the one determined using single plane super resolution image is now shown on Figure 3B-G. We have also added a new Figure 3—figure supplement 1 with further images and a quantification of the differences between polarities determined with the two approaches. The quantification shows that more than 80% of PIN1 cellular polarities match almost perfectly (deviation less than 30°) when comparing the two approaches (Figure 3—figure supplement 1 and corresponding text). We have also reorganized the supplementary material to include a Supplementary Methods 3 that is dedicated specifically to the analysis of PIN1 polarity and provide extensive details on the methods we have used. Note that the approach we describe in the manuscript allows to obtain not only quantitative information on the direction of the PIN1-mediated transport but also on the intensity of the protein on each wall of the meristem, which is a likely determinant parameter of the transport intensity.

Also we have removed the sentence pointed out by the reviewer, as it was unnecessarily vague and, in doing so, we keep the focus on the parameters that we have analyzed in order to keep the claim of the manuscript in line with the data we present. For the same reason, we also insist on the fact that we have aimed at providing a quantitative tissue-scale vision of the organization of PIN1 polarities and we have been careful to avoid premature conclusions.

For the localization of the auxin sources, we provide details in our response to the point 3 of the reviewer below.

2) Even if one could measure all necessary parameters at the desired resolution, to make a point that the new findings lead to a consistent new picture, extremely detailed mathematical modelling would have to accompany each proposed novel mechanism.

It is fair to say that the consistency of the new picture we obtain would ultimately require mathematical modeling. While we are currently working on such models, adding those to the manuscript effectively goes beyond what can be reasonably included in the manuscript. Again our goal on the part concerning auxin transport and biosynthesis is to provide the community with a high-resolution tissue scale spatio-temporal reference dataset with information on auxin, distribution, its transport and its biosynthesis.

3) As even the best imaging technology merely establishes reliable correlation, there has always been criticism when studies were entirely based on imaging. There is no genetics and very few application experiments in this study to critically test the modelled hypotheses.

To address this concern, we have now added genetic data to support two conclusions of the manuscript. This first concerns the role of YUC and biosynthesis on SAM function. Genetic data have been previously published (Shi et al., 2018; Pinon et al., 2013; Cheng 2006); our expression data motivated the addition in the revised manuscript of an analysis of the yuc1yuc4 mutants. yuc1yuc4 double mutants show organ positioning and morphological defects as shown on Figure 3L and Figure 3—figure supplement 3O-U. In addition to published work, this support the idea that auxin biosynthesis from YUC1 and YUC4 in flowers is important for phyllotaxis. This is explained in the text in the following way: “In addition, yuc1yuc4 loss-offunction mutants show severe defects in SAM organ positioning and size (Shi et al., 2018; Pinon et al., 2013) (Figure 3L and Figure 3—figure supplement 3O-U). Taken with the organization of PIN1 polarities, these results suggest that P3-P5 are auxin production centers for the SAM that regulate phyllotaxis.”

The second conclusion concerns time integration of auxin information in the regulation of transcription during organogenesis. We have now analyzed the effect of auxin treatments on the induction of transcription using different concentrations and durations of treatment. The results are included in Figure 5I that shows the effect on DR5 expression of treatments with 0.2mM, 1mM and 5 mM auxin for different durations. This allows us now to demonstrate that both the time of exposure and the concentration of auxin control the activation of transcription monitored by DR5, thus supporting the conclusion that cells integrate the concentration of auxin over time and not simply time. In addition, we have analyzed transcription activation by auxin (using DR5) in pinshaped inflorescence of the pinoid mutant. Organ initiation is inhibited in pinoid meristems and our results indicate that there is no pre-established auxin distribution. We also show that all cells at the periphery of the meristem can respond similarly to an exogenous auxin treatment, which provides further argument that there are no intrinsic differences in the capacity of the cells to respond to auxin. We have included these new results in Figure 5J-M and Figure 5—figure supplement 2A-C. We explain the effect of auxin treatments in wild-type and pinoid meristems in the text.

Other reviewers then suggested that genetics could be used to show an involvement of the auxin signaling pathway in time integration of the auxin signal. There are few auxin signaling proteins for which a clear function has been established in the SAM. The one with the clearest implication is ARF5/MP but the mutant was shown to be insensitive to auxin in the meristem (Reinhardt et al., 2003) and thus cannot be used here. ARF3/ETTIN and ARF4 are the two other ARFs that show high expression during organ initiation (Vernoux et al., 2011). A key role for ARF3/ETTIN in organ initiation at the SAM (alongside ARF4) has recently been demonstrated by a publication from the laboratory of Doris Wagner. We have included in the revision data demonstrating that the ARF3/ETTIN gene is necessary for the integration of the auxin signal. These new data are included in Figure 5 (Figure 5N-Q) and a new Figure 5—figure supplement 2 (Figure 5—figure supplement 2D,E) and described in the text in the following way:

“The Auxin Response Factor (ARF) ETTIN (ETT/ARF3) plays an important role in promoting organogenesis in the SAM (Wu et al., 2015; Chung et al., 2019). We found that in a loss-of function ett3 mutant the expression of DR5 was restricted to only 2-3 cells at sites of organogenesis, and that a 300’ 1mM auxin treatment did not induce DR5 in the SAM (Figure 5N-Q and Figure 5—figure supplement 2D-E).”

As the functional link between auxin signaling and histone acetylation is well established,, we now show the results of TSA treatments on DR5 expression in a supplementary figure (Figure 5—figure supplement 1) and discuss them in the following way after the description of the result with ARF3:

“Auxin signaling and ARF3 have been shown to act by modifying acetylation of histones (Wu et al., 2015; Chung et al., 2019 and Long et al., 2006). Pharmacological inhibition of histone deacetylases (HDACs) alone was able to trigger concomitant activation of DR5 at P0 and P-1 sites in the SAM (Figure 5—figure supplement 1I-K). These results suggest that auxin signal integration depends on a functional ARF-dependent auxin nuclear pathway.”

Here, the combination of auxin treatments and genetics allows us to provide a deeper understanding of how spatial and temporal auxin information is used to trigger transcription of target genes during organ initiation.

4) While being detailed in experimental documentation, at the conceptual level, the paper is rather poorly written. The following sub-points may (incompletely) illustrate, how this manuscript text would fall short of the expectations, for example in a developmental genetics graduate course.

– What are the authors distinguishable definitions for “morphogen” vs. “morphogenetic regulators“?

– The selection of references is correspondingly arbitrary. That substances “trigger” certain cell behaviors in some of the references is certainly insufficient to put them in the “morphogen/morphogenetic” category. Likewise, all kinds of gradients are abundant in biology; do the authors claim that distinct concentrations of auxin correspond to distinct cell fates or not?

Auxin has not been formally established as a morphogen according to the criteria used in animal systems. Auxin is clearly a regulator of morphogenesis and we used the term “morphogenetic regulator” in the manuscript for this only reason. But we agree that this introduced an unnecessary confusion. Addressing whether auxin indeed acts as a morphogen goes beyond the scope of the manuscript, although we have established in this work some key tools to do so. We have thus removed the term “morphogenetic regulator” from the revised manuscript and rather focus on the concept of positional information. Corrections have been introduced throughout the manuscript and we believe that we now provide a solid text at the conceptual level.

– Because of the pivotal role the epidermis had in earlier SAM modelling, it is surprising that the main text does not address this issue. To mention in the technical sections that an experiment “focuses” on the epidermis, does not answer what conceptual roles different tissue layers have.

The pivotal role of the epidermis is also established genetically, with the demonstration that the pin1 mutant can be complemented through expression of PIN1 only in the epidermis (Kierzkowski et al., 2013). We have clarified this in the main text: “We used the spatial distribution of 1-DII-VENUS/TagBFP as a proxy for auxin distribution, hereafter named “auxin” (Figure 1C) and focused on the epidermal cell layer (L1) where organ initiation takes place (Jonsson et al., 2006; Kierzkowski et al., 2013; Smith et al., 2006; Reinhardt et al., 2003).”

Limited genetic and functional data make it difficult to conceptualize the function of auxin in the different tissue layers and the respective contribution of these layers. As the manuscript’s focus is on linking auxin information to transcriptional activation, which also occurs first in the L1 layer, analyzing the L1 layer is logical. Addressing the role of the different tissue layers would require to analyze specifically this question, which goes beyond the scope of this manuscript.

– Some of the co-authors have implicated auxin, along with cytokinin in shaping the SAM stem cell zone (Leibfried et al., 2005, Zhao et al., 2010). How are those multiple roles of the same signal substance(s) are integrated conceptually?

We agree that this is an interesting question. There is a paragraph in the discussion (that was already in the previous version of the manuscript) that discusses how the dynamics of cell differentiation when cells exit the stem cell niche might be regulated through opposite effects of auxin and WUS on the histone acetylation status of auxin signaling and target genes. In this paragraph, we refer the reader to a recent paper some of the co-authors published in Nature Com (Ma et al., 2019). This paper further demonstrates the need of a low level of auxin signaling for stem cell activity. It is possible that different levels of auxin signaling activity (and possibly other regulations) lead to activation of different set of genes in the PZ and CZ, explaining different roles. However a significant amount of work is still necessary to understand this. As our work focuses primarily on the action of auxin in organogenesis, we have not added extra discussion on this point in the revised manuscript.

Second round of Review

Reviewer 1:

In this revision many of my previous concerns were addressed, notably controls about the timing vs. total amount of auxin application and a de-emphasis of the HDAC results.

This new version keeps it strengths in quantitative imaging and reveals a number of fascinating phenomenon: namely the high degree to which auxin (as read out by DIIVenus) accumulates in stereotyped patterns in the SAM and primordia and the fact that these patterns are not merely a consequence of cell and tissue growth. The model of time integration for spatial outputs, while not completely nailed down, is interesting and supported well enough to be a reasonable hypothesis.

We note that the reviewer still finds that we report fascinating phenomenons.

There are still some sections where I do not follow the logic and/or feel that the conclusion is not justified by the data. Specifically the ETT/HDAC section is not particularly well explained or convincing. Inhibition of HDACs by broad chemical is likely to be different than the behavior of ETT toward target histones. To make the case that ETTs works through HDAC I would expect to see WT, ET + HDAC inhibted and ett mutant meristems probed for the specific histone acetylation and a correlation between the affected cells/ regions and auxin behaviors. Line 258-9, is there an alternative to the idea that auxin signal integration requires an auxin dependent nuclear pathway? And how does inhibiting HDACs show this?

DR5 and DII-Venus may be the most commonly used markers of auxin levels and transcriptional response, but they are still proxies, and in the case of DR5, almost certainly inaccurate in the specifics. I think the authors should refer to DII-VENUS the same way they do DR5—by name rather than “auxin” or “auxin input”. In my view, the degree to which the data presented in Figure 4 represents TRUE auxin input and output is still up for debate, but I think, given the way the community has used these reporters, the evidence of them not correlating in this assay is useful.

We feel that most of these comments can be easily addressed though adjustment of the text. Although this is stated in the text, we might have not been fully clear on the fact that the demonstration of a role for HDAC in regulation of genes by ETT has already been well established by the team of Doris Wagner in Chung et al., 2019. We thus feel that it is not needed to “make the case that ETT works through HDAC”, in accordance with the fact that the reviewer appreciates that we “de-emphasized the HDAC results”.

Reviewer 4:

The manuscript by Galvan-Ampudia et al. is a resubmission of a manuscript I had previously reviewed. In that review, I suggested that the Authors had provided convincing evidence that in the epidermis of the inflorescence apex of Arabidopsis cells with peak expression of a proxy of auxin concentration do not move apart because of intercalating cell growth or proliferation; rather, it's the maxima themselves that move across cells over time. My biggest concern with the previous version of this manuscript was the lack of biological relevance for the claim — derived from the observations above — that it's the exposure of cells to auxin for a certain amount of time that defines how those cells will respond to auxin.

I believe Figure 5A–I convincingly shows that treatment with 0.2 mM IAA for 300 minutes leads to – at least – the same level of DR5 expression that treatment with 1 mM IAA for 120 minutes does. However, the former treatment only led to 50% of plants with defective phyllotaxis, whereas the latter treatment gave rise to nearly 90% of plants with defective phyllotaxis (Figure 5—figure supplement 1D). Therefore, shorter treatment with higher auxin concentration led to nearly twice as many abnormal plants than longer treatment with higher auxin concentration, suggesting that the amount of time cells are exposed to auxin is perhaps not as biologically relevant as the Authors wish to claim; this reduces my enthusiasm for this manuscript.

I was suprised to see that an important experiment such as the one reported in Figure 5—figure supplement 1D — one that is meant to provide evidence of the biological relevance of the Authors' observations — had been relegated to a single panel in the supplemental material. I was also suprised to find that both the number of plants on which that experiment had been performed and the statistical analysis of those data were missing. And I was surprised to see that the Authors grouped all plant responses to the treatments into a single, uninformative phenotype class vaguely defined as "defective phyllotaxis". This latter deficiency may have precluded the identification of components of plant response to the treatment that are specfically affected by the time cells are exposed to auxin. For example, though the fraction of plants showing "phyllotaxis defects" is lower when plants are treated with 0.2 mM IAA for 300 minutes than when they are treated with 1 mM IAA for 120 minutes, it is possible that the fraction of plants with a specific defect is instead higher. This could suggest what biological parameter, if any, depends on temporal integration of auxin signals; most important, it would also suggest that such temporal integration does matter for plant growth and development.

Adjusting the figure can easily be done if needed. The reviewer is right that the effects of auxin treatments could raise questions when different concentrations are compared. However, we would like to point that a duration-dependent effect is seen for the different concentrations and that the reviewer seems to have overlooked our analysis of the ett mutant that provide evidence at the SAM level of the effect of an alteration of time integration of the auxin signal. For this reason, we disagree with the final conclusion of the reviewer. But it could be that the way we did our analysis of phyllotaxis are the reason for the inconsistencies. We have some more detailed information on this experiment that we are currently analyzing. We have also plants almost ready to repeat this experiment in depth if needed.

On a minor note, the response to the reviewers' comments I had access to mentioned the use of the MP::MPDelta:GR and MP::MPDelta:EAR:GR lines generated by the Jiao lab, but I could not find any sign of them in the resubmitted manuscript.

We had simply problems to have this material working in our hands and the analysis of the ett mutant gave us most of the answers we needed.

Reviewer 5:

The submitted study presents novel reporter constructs and groundbreaking imaging work that sheds new light on the distribution of auxin during organ formation in the shoot meristem. In particular the authors show precise temporal dynamics and the impressive reproducibility of auxin distribution patterns, plus they chart those relative to the movements of the participating cells. Not least, they relate those distributions to patterns of transcriptional responses to auxin, which could reflect a cellular memory of auxin exposure and the system is interrogated by external auxin applications of different concentrations and durations. Finally, the manuscript provides data on the dynamics of PIN1 polarities and maps expression sites of YUCCA genes. The strength of the manuscript clearly lies in the novel and superior precision of auxin distribution recordings. Moreover, the authors show highly reproducible dynamic (temporal) patterns with interesting relationships to the dynamics of cell proliferation in the same areas. These data provide solid experimental ground to anchor future and to test current mathematical models that link auxin distribution patterns and cell behaviour in shoot meristems. There is also some progress on quantification of the “hysteresis” effect, although the manipulations are relatively crude. Together, these findings constitute a solid and important contribution and a more relaxed density of presentation would allow for the integration of Figure 1—figure supplement 2 into the body of the manuscript, as it is too central for being left out.

We note that the reviewer praised the solidity of the findings. We would like however to point that we do not understand why our manipulations are “crude”. We did extensive quantifications of transcriptional activation down to a cellular level in treatments with 3 auxin concentrations and 2-3 durations of treatments for each concentration. We report precisely how transcription changes relatively in space, with adequate statistical processing to support our conclusions. We thus believe that we have analyzed in depth the “hysteresis effect”. Concerning Figure 1—figure supplement 2, it could certainly be a main figure if needed.

It is not really clear why those solid data should to be diluted by weaker data just to generate a comprehensive model, which would be as weak as its weakest support. The PIN polarity data are not fundamentally superior to other published records and their functional relevance as well as those of YUCCA expression data is not experimentally challenged. There are only two (as it appears, arbitrarily chosen) genotypes introduced, although there would be numerous genotypes available, including (but by far not limited to) those enabling local targeted manipulation of auxin distribution and response. Such genetic tests could challenge and provide functional support for any model.

Here, the reviewer might have overlooked the fact that we have used super-resolution microscopy to co-visualize PIN1 and PI staining in order to validate our quantitative analysis of PIN1 polarities using confocal microscopy. This is superior to analysis using visual inspection (as found in the literature) that do not provide quantitative data. The reviewer might also have overlooked the fact that we have experimentally challenged YUCCA expression using mutants (in complement of published studies) and that these mutants were not arbitrarily chosen but chosen on expression as we provided an analysis of the expression of all YUCCAs in the SAM. Our results could effectively be further tested using genetics but we believe that it is out of the scope of the manuscript.

https://doi.org/10.7554/eLife.55832.sa2

Article and author information

Author details

  1. Carlos S Galvan-Ampudia

    Laboratoire Reproduction et Développement des Plantes, Univ Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Inria, Lyon, France
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing
    Contributed equally with
    Guillaume Cerutti
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4779-3568
  2. Guillaume Cerutti

    Laboratoire Reproduction et Développement des Plantes, Univ Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Inria, Lyon, France
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing
    Contributed equally with
    Carlos S Galvan-Ampudia
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2003-9335
  3. Jonathan Legrand

    Laboratoire Reproduction et Développement des Plantes, Univ Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Inria, Lyon, France
    Contribution
    Resources, Software, Formal analysis, Validation, Investigation, Methodology
    Competing interests
    No competing interests declared
  4. Géraldine Brunoud

    Laboratoire Reproduction et Développement des Plantes, Univ Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Inria, Lyon, France
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
  5. Raquel Martin-Arevalillo

    Laboratoire Reproduction et Développement des Plantes, Univ Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Inria, Lyon, France
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
  6. Romain Azais

    Laboratoire Reproduction et Développement des Plantes, Univ Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Inria, Lyon, France
    Contribution
    Formal analysis, Validation
    Competing interests
    No competing interests declared
  7. Vincent Bayle

    Laboratoire Reproduction et Développement des Plantes, Univ Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Inria, Lyon, France
    Contribution
    Formal analysis, Investigation
    Competing interests
    No competing interests declared
  8. Steven Moussu

    Laboratoire Reproduction et Développement des Plantes, Univ Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Inria, Lyon, France
    Present address
    Department of Plant Molecular Biology, University of Lausanne, Lausanne, Switzerland
    Contribution
    Methodology
    Competing interests
    No competing interests declared
  9. Christian Wenzl

    Department of Stem Cell Biology, Centre for Organismal Studies, Heidelberg University, Heidelberg, Germany
    Contribution
    Resources
    Competing interests
    No competing interests declared
  10. Yvon Jaillais

    Laboratoire Reproduction et Développement des Plantes, Univ Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Inria, Lyon, France
    Contribution
    Conceptualization, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4923-883X
  11. Jan U Lohmann

    Department of Stem Cell Biology, Centre for Organismal Studies, Heidelberg University, Heidelberg, Germany
    Contribution
    Resources, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3667-187X
  12. Christophe Godin

    Laboratoire Reproduction et Développement des Plantes, Univ Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Inria, Lyon, France
    Contribution
    Conceptualization, Supervision, Funding acquisition, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    christophe.godin@ens-lyon.fr
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1202-8460
  13. Teva Vernoux

    Laboratoire Reproduction et Développement des Plantes, Univ Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Inria, Lyon, France
    Contribution
    Conceptualization, Supervision, Funding acquisition, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    teva.vernoux@ens-lyon.fr
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8257-4088

Funding

Agence Nationale de la Recherche (ANR-12-BSV6-005)

  • Teva Vernoux

Deutsche Forschungsgemeinschaft (FOR2581)

  • Christian Wenzl
  • Jan U Lohmann

Human Frontier Science Program (RPG0054-2013)

  • Christophe Godin
  • Teva Vernoux

Agence Nationale de la Recherche (ANR-18-CE12-0014-02)

  • Teva Vernoux

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

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 and mutant lines; Sophie Ribes for her contribution to the nuclei image processing pipeline; Gwyneth Ingram, Olivier Hamant, Arezki Boudaoud and Jan Traas for feedback on the manuscript. We acknowledge the contribution of SFR Biosciences (UMS3444/CNRS, US8/Inserm, ENS de Lyon, UCBL) facilities PLATIM, for assistance with microscopy; of the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for providing us access to their computing resources. This work was supported by Human Frontier Science Program organization (HFSP) Grant RPG0054-2013 to TV and CG, ANR-12-BSV6-0005 grant (AuxiFlo) and ANR-18-CE12-0014-02 (ChromAuxi) to TV and DFG FOR2581 to JUL.

Senior Editor

  1. Christian S Hardtke, University of Lausanne, Switzerland

Reviewing Editor

  1. Jürgen Kleine-Vehn, University of Natural Resources and Life Sciences, Austria

Publication history

  1. Received: February 7, 2020
  2. Accepted: March 29, 2020
  3. Version of Record published: May 7, 2020 (version 1)

Copyright

© 2020, Galvan-Ampudia et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,810
    Page views
  • 415
    Downloads
  • 2
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Developmental Biology
    Elisa Nerli et al.
    Research Article Updated

    During brain development, progenitor cells need to balanceproliferation and differentiation in order to generate different neurons in the correct numbers and proportions. Currently, the patterns of multipotent progenitor divisions that lead to neurogenic entry and the factors that regulate them are not fully understood. We here use the zebrafish retina to address this gap, exploiting its suitability for quantitative live-imaging. We show that early neurogenic progenitors arise from asymmetric divisions. Notch regulates this asymmetry, as when inhibited, symmetric divisions producing two neurogenic progenitors occur. Surprisingly however, Notch does not act through an apicobasal activity gradient as previously suggested, but through asymmetric inheritance of Sara-positive endosomes. Further, the resulting neurogenic progenitors show cell biological features different from multipotent progenitors, raising the possibility that an intermediate progenitor state exists in the retina. Our study thus reveals new insights into the regulation of proliferative and differentiative events during central nervous system development.

    1. Cell Biology
    2. Developmental Biology
    Jamie C Little et al.
    Research Article Updated

    Extracellular Hedgehog (Hh) proteins induce transcriptional changes in target cells by inhibiting the proteolytic processing of full-length Drosophila Ci or mammalian Gli proteins to nuclear transcriptional repressors and by activating the full-length Ci or Gli proteins. We used Ci variants expressed at physiological levels to investigate the contributions of these mechanisms to dose-dependent Hh signaling in Drosophila wing imaginal discs. Ci variants that cannot be processed supported a normal pattern of graded target gene activation and the development of adults with normal wing morphology, when supplemented by constitutive Ci repressor, showing that Hh can signal normally in the absence of regulated processing. The processing-resistant Ci variants were also significantly activated in the absence of Hh by elimination of Cos2, likely acting through binding the CORD domain of Ci, or PKA, revealing separate inhibitory roles of these two components in addition to their well-established roles in promoting Ci processing.