# Robust branch patterning in moss shoots via symplasmic auxin diffusion

Jeanne Abitbol-Spangaro, Gladys Cloarec, Arthur Muller, Stéphanie Hallet, Claire Boulogne, Cynthia Gillet, Vojtěch Schmidt, Petre Dobrev, Roman Skokan, Valentin Couvreur, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Jeanne Abitbol-Spangaro, Gladys Cloarec, Arthur Muller, Stéphanie Hallet, Claire Boulogne, et al. Robust branch patterning in moss shoots via symplasmic auxin diffusion. Current Biology, In press, 10.1016/j.cub.2025.09.031. hal-05306843

HAL Id: hal-05306843

https://hal.science/hal-05306843v1

Submitted on 9 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Robust branch patterning in moss shoots
via symplasmic auxin diffusion

Jeanne Abitbol-Spangaro<fnref n="1" />, Gladys Cloarec<fnref n="1" />,<fnref n="6" />, Arthur Muller<fnref n="1" />,<fnref n="6" />, Stéphanie Hallet<fnref n="1" />,<fnref n="6" />, Claire Boulogne<fnref n="2" />, Cynthia Gillet<fnref n="2" />, Vojtěch Schmidt<fnref n="3" />,<fnref n="4" />, Petre I. Dobrev<fnref n="3" />, Roman Skokan<fnref n="3" />, Valentin Couvreur<fnref n="5" />, Jeroen de Keijzer<fnref n="6" />, Christophe Godin<fnref n="1" />,<fnref n="7" />, and Yoan Coudert<fnref n="1" />,<fnref n="7" />,<fnref n="8" />

## 1 Summary

Uncovering the mechanisms by which developmental patterns consistently arise within species is a major goal of biology. In plants, branching patterns are key determinants of morphological diversity. Similar lateral branching modes convergently evolved in the leafy shoot of vascular plants and bryophytes, driving their independent architectural diversification. While auxin-dependent branch inhibition is shared between both lineages, long-range polar auxin transport plays a central role in branching control in vascular plants, but in bryophytes like the moss *Physcomitrium patens*, auxin might diffuse through plasmodesmata to regulate branch distribution. However, whether symplasmic auxin diffusion is a biophysically realistic mechanism sufficient to explain observed branch distribution patterns has yet to be assessed. Here, we address this fundamental problem by developing a physics-based, three-dimensional computational model of symplasmic auxin diffusion in the moss shoot, integrating molecular, cell and tissue scales. Each step of model design was guided by geometry measurements and biological experiments. Our integrative approach demonstrates that branching control based solely on symplasmic diffusion for intercellular auxin movement can account for the observed branching patterns at the whole-shoot level. It also provides mechanistic interpretations of the changes in branch distribution caused by genetic perturbations affecting callose-dependent symplasmic permeability, as well as the unexpected increase in branch spacing robustness during shoot development. Altogether, our findings reveal that branching patterns arising from an auxin diffusion-based regulatory mechanism exhibit a specific developmental signature, not reported in vascular plants, but well exemplified in the moss *Physcomitrium*.

**Keywords :** plasmodesmata, moss, *Physcomitrium*, branching, auxin, diffusion

## 19 Introduction

Developmental patterns are defining features of multicellular organisms. A major goal of biology is to understand how these patterns consistently emerge within species, by identifying the molecular regulatory signals involved and

<footnote n="1">¹Laboratoire Reproduction et Développement des Plantes, ENS de Lyon, CNRS, Université Claude Bernard Lyon 1, INRAE, INRIA, 69007 Lyon, France</footnote>

<footnote n="2">²Electron Microscopy Facility, Imagerie-Gif, Institute for Integrative Biology of the Cell (I2BC), CEA, CNRS, Univ. Paris-Sud, Université Paris-Saclay, 91198 Gif-sur-Yvette cedex, France</footnote>

<footnote n="3">³Institute of Experimental Botany of the Czech Academy of Sciences, Rozvojová 263, 165 00 Prague 6, Czechia</footnote>

<footnote n="4">⁴Department of Experimental Plant Biology, Charles University, Viničná 5, 128 00 Prague 2, Czechia</footnote>

<footnote n="5">⁵Université Catholique de Louvain, Earth and Life Institute - Agronomy, 1348 Louvain-la-Neuve, Belgium</footnote>

<footnote n="6">⁶Laboratory of Cell and Developmental Biology, Wageningen University and Research, Droevendaalsesteeg 1, 6708 PB Wageningen, The Netherlands</footnote>

<footnote n="6">⁶Denotes equal contribution</footnote>

<footnote n="7">⁷Correspondence : christophe.godin@inria.fr, yoan.coudert@cnrs.fr</footnote>

<footnote n="8">⁸Lead contact</footnote>
the mechanisms by which they are integrated across cells and tissues¹. In plants, branch distribution is a key developmental pattern contributing to their morphological diversity²,³. Over the past century, research in vascular plant models has established the phytohormone auxin as a major regulator of branching patterns⁴⁻⁷. In particular, the formation of auxin minima in the axils of developing leaves is required to promote the initiation of branch meristems, *i.e.*, the founder cells that may later develop into lateral branches. Additionally, auxin biosynthesis at the shoot apical meristem is also crucial for maintaining apical dominance, the process by which a growing apex inhibits the development of distant buds or branches⁸⁻¹⁴.

Recently, studies using the moss *Physcomitrium patens* as a model have revealed that auxin’s inhibitory role in determining branch position is shared across land plants. Specifically, branching in the moss gametophore (*i.e.*, the leafy shoot) likely occurs in a single step where new meristems form only in the axils of certain developed phyllids (*i.e.*, the leaf-like organs) and directly grow into branches, in response to specific signals¹⁵. Auxin is produced at the base and in the apices of both gametophores and their lateral branches. *Physcomitrium* mutants lacking *TRYPTOPHAN AMINOTRANSFERASE (TAR)* or *YUCCA FLAVIN MONOOXYGENASE-LIKE (YUC)* auxin biosynthesis genes exhibit reduced apical dominance and increased gametophore branching, suggesting that auxin inhibitory action requires its movement away from sources¹³. However, in contrast with vascular plants, where PIN-FORMED (PIN)-mediated basipetal polar auxin transport is a central regulatory mechanism of shoot branching, gametophores lack long-range polar auxin transport and PIN auxin efflux transporters are not essential branching regulators in *Physcomitrium*¹⁰,¹⁵⁻¹⁷. Supporting this view, auxin-mediated inhibition is sufficient to explain *Physcomitrium* branching patterns (which are characterized by apical, inter-branch and basal zones devoid of branches, also called branch inhibition zones) only if auxin circulates in an apolar manner from its sources¹⁵. These results have led to the hypothesis that intercellular auxin movement in gametophore stems may occur by symplasmic diffusion, *i.e.*, through plasmodesmata.

Although this pathway has yet to be proven sufficient for branching regulation in moss, genetic studies in other systems have demonstrated that auxin can effectively diffuse between cells through plasmodesmata¹⁸. Symplasmic permeability is modulated by the deposition and turn-over of callose, a β-1,3-glucan, at the neck of plasmodesmata¹⁹. This process is regulated by evolutionarily conserved CALLOSE SYNTHASE (CALS) and GLYCOSYL HYDROLASE-LIKE FAMILY 17 (GHL17) enzymes that promote callose synthesis and degradation, respectively²⁰. *Arabidopsis* plants with altered CALS or GHL17 function exhibit modified auxin distribution patterns in various developmental contexts, indicating that callose-dependent plasmodesmal gating controls auxin movement²¹⁻²⁴. However, both plasmodesmata-mediated auxin diffusion and active auxin transport pathways are generally required together to explain the observed auxin distribution across tissues²⁵,²⁶, leaving open the question of whether symplasmic diffusion alone could be sufficient to account for auxin-based branch patterning.

To further evaluate the plausibility of this hypothesis in *Physcomitrium*, it is necessary to verify that the characteristic diffusion length through plasmodesmata matches the observed size range of the macroscopic inhibition mechanism. This is of critical importance since molecular diffusion is a slow mechanism, operating at a local scale, which alone is often insufficient to generate robust developmental patterns (*i.e.*, consistent phenotypic outcomes despite variable inputs)²⁷. For example, diffusion is sufficient to establish the Fibroblast growth factor 8 (Fgf8) gradient controlling zebrafish embryo development because the required characteristic length is about 200µm²⁸. Nevertheless, in biological systems requiring morphogenetic control over larger scales, additional transport or relay mechanisms are necessary²⁹. In plants, the characteristic length scale of auxin diffusion in an idealized, continuous and one-dimensional tissue can be estimated to $\lambda = \sqrt{D/\delta} \approx 1\text{mm}$ (with auxin diffusivity in the cytoplasm, $D = 300\text{µm s}^{-130}$; and auxin decay rate in *Arabidopsis*, $\delta = 1\text{h}^{-131}$). However, auxin diffusivity in real, three-dimensional (3D) plant tissues is likely less due to the physical properties of the symplasmic continuum, which are notably determined by plasmodesmata density and permeability. Additionally, auxin diffusivity may be influenced by the overall organ geometry, which can reflect or absorb auxin’s random movements at organ boundaries. Effective
auxin diffusion in complex multicellular structures is therefore a composite variable, arising from the interplay of various physiological and structural components across multiple scales.

The aim of this work is to evaluate whether the combination of these components in the context of the 3D moss gametophore can produce a realistic diffusive mechanism explaining branching patterns. To address this, key geometric and biological properties have been estimated at multiple scales, including plasmodesmata density and permeability, auxin physiology, and cell, tissue and organ-level 3D geometry. These quantified constituents have been integrated into a 3D computational model to simulate auxin diffusion within the gametophore's symplasmic network. The plausibility of the model has then been tested based on its ability to recapitulate observed branch inhibition zones in *Physcomitrium*, predict responses to genetic perturbations affecting symplasmic permeability, and reveal remarkable emergent properties of this original branching mechanism.

## Results

**Modeling symplasmic auxin diffusion in the moss gametophore.** To construct our model, we first created a simplified 3D representation of the gametophore stem's geometry (Figures 1A and 1B). We began with an idealized cross-section comprising cells generated through a Voronoi tesselation of a disk and extruded this cross-section over several cell layers of varying heights. The resulting 3D cellular structure consisted of polyhedral cells (Figure 1C).

To model auxin movement in this virtual tissue, we assumed that there was no carrier-mediated active auxin transport in the gametophore. In addition, we considered that auxin diffusion was primarily limited at cell interfaces by the surface area of plasmodesmata openings available for cell-to-cell exchange, making intracellular auxin diffusion effectively instantaneous by comparison. At cell interfaces, auxin was assumed to diffuse across the cytoplasmic sleeve of plasmodesmata, thus with typical cytosolic diffusivity, but through a restricted surface area shaped by the local geometrical and topological properties of the symplasmic continuum. Cell face permeability was therefore defined as the product of the plasmodesmata permeability $P$ and the plasmodesmata number. Plasmodesmata permeability $P$ was calculated as the product of auxin's molecular diffusivity in the cytosol $D$, the plasmodesmata open-cross section $\omega$ and the inverse of the cell face thickness $l$, primarily corresponding to the cell wall thickness. The plasmodesmata number was defined as the product of the plasmodesmata density $\rho(i, j)$ at the interface between cells $i$ and $j$ and the interface area $A_{i,j}$ (Figure 1D). Diffusion fluxes were considered proportional to auxin concentration differences between adjacent cells, according to Fick's first law. We also assumed that auxin decay occurred in all cells at a constant rate $\delta$.

Based on these assumptions, the temporal variation of auxin concentration in each cell was expressed as the sum of degradation and diffusion processes, using a set of coupled ordinary differential equations applied to the graph representing the tissue (Figure 1D), such that for $i = 1...N$,

$$ \frac{da_i(t)}{dt} = \frac{1}{V_i} \sum_{j \in N_i} P A_{i,j} \rho(i, j) (a_j(t) - a_i(t)) - \delta a_i \quad (1) $$

with $V_i$ the cell volume. Additionally, auxin concentration was fixed at a constant value $\alpha$ in a set of source cells $S$, such that $a_s(t) = \alpha$, $\forall s \in S$, $\forall t \ge 0$, and non-source cells initially contained no auxin, i.e., $a_i(0) = 0$, $\forall i \notin S$.

**Calibration and estimation of model parameters.** We assumed that parameters $D$, $l$, $\omega$ and $\delta$ were constant in the model, and $D$, $l$ and $\omega$ values were estimated from published studies (Data S1A). However, the rate of auxin decay $\delta$ was unknown in *Physcomitrium*. To estimate $\delta$ value, gametophores were supplied with $[^{13}C_6]$IAA, a stable isotope-labeled form of auxin, and the incorporation of $[^{13}C_6]$ into known auxin catabolites was tracked over time (Figure S1). We found that IAA was predominantly converted into 2-oxindole-3-acetic acid (oxIAA), indicating that oxIAA is the dominant auxin catabolite in *Physcomitrium*, consistent with previous findings in other moss
species³². We then calculated the rate of IAA-to-oxIAA conversion as an approximation for the auxin decay rate, yielding an estimated value of 3h⁻¹, which represents an upper bound of the decay rate given our experimental setup.

Another critical parameter of the model was the plasmodesmata density ρ(i, j). To determine plasmodesmata densities across the entire gametophore, whole gametophores were fixed, transversely sectioned from the apex towards the base, and relevant sections were selected for electron microscopy imaging (Figures 2A, S2 and S3). Cell walls were imaged at high resolution to count individual cell wall-traversing cytoplasmic strands, hereafter referred to as plasmodesmata, and map their frequencies across entire stem sections. We found that young and undifferentiated tissues located at the gametophore apex were rich in plasmodesmata, with an average density of 35 µm⁻² (Figures 2B and 2E). Away from the apex, in regions with differentiated tissues as suggested by cell vacuolation, plasmodesmata densities decreased to an average value of 1 µm⁻² (Figures S3C-S3E and S3I). Nevertheless, in these regions, plasmodesmata densities were comparatively much higher in branch meristem cells than in surrounding tissues (Figures 2C, 2E, 2F and S3I). Cells at the base of well-grown gametophore branches were well connected with an average plasmodesmata density of 2.4 µm⁻² (Figures S3G-S3I). To complement these observations, a gametophore was longitudinally sectioned through its central axis, which confirmed that plasmodesmata densities were higher around the apex (30 µm⁻²) (Figure 2D). This also revealed more elevated densities on transverse than longitudinal walls (7.8 versus 4.2 µm⁻² on average, respectively), suggesting that plasmodesmata on longitudinal walls form near the apex and their densities diminish during growth, likely resulting from a dilution process caused by wall elongation (Figures 2G, 2H, S3J and S3K). Thus, at the gametophore level, plasmodesmata distribution is characterized by density gradients that peak in meristematic regions and decrease as cells differentiate.

Gametophore geometry was also a key aspect of our model. To achieve realistic geometric dimensions in our 3D tissue representation, we measured the average transverse cell wall length (13 µm) and cell area (400 µm²) from transverse gametophore sections and used these measurements to define cell sizes in the idealized cross-section (Figures 3A and S4A-S4C). Then, to define the thickness of cell layers based on their distance from the apex, we quantified cell length along the apical-basal gametophore axis, which varied from 10 to 200 µm with a median of 80 µm in the mature zone (Figures 3B and S4D-S4F). The length of phytomers in the mature part of the gametophores was estimated at 236 µm on average, suggesting that mature phytomers consisted of an average of three cells (Figures S4G and S4H). This was consistent with direct microscopic observations of the gametophore surface (Figure S4I). In addition, based on previous observations of phyllotaxis in *P. patens*³³, we assumed that phyllids were separated by an average angle of 137.5° in the azimuthal direction on the gametophore (Figure S4J).

Finally, we modeled the observed plasmodesmata distribution within our cellular framework. This was achieved using a Hill function h dependent on the distance to the auxin source, assuming that plasmodesmata dilution was proportional to the face area for longitudinal cell faces. Specifically, $d_{ij}$ represented the distance in µm from cell face (i, j) to the closest auxin source, such that:

$$ \rho(d_{ij}) = \begin{cases} h(d_{ij}) & \text{for transverse cell faces} \\ h(d_{ij})A_{trans}/A_{ij} & \text{for longitudinal cell faces} \end{cases} \qquad (2) $$

where $A_{trans}$ was the average area of transverse cell faces. The Hill function h was optimized to capture the plasmodesmata densities observed in the different gametophore sections (Figures 3C-3E and S5, see STAR Methods). This approach also allowed us to account for the absence of certain features of the biological tissues in the model, such as staggered cells and a conical apex. Since $d_{ij}$ was independent of tissue dimensionality, the plasmodesmata distributions estimated in 2D tissues could be seamlessly applied to the 3D cellular framework.

Collectively, these data allowed us to construct a 3D model incorporating realistic gametophore tissue and cell geometry, and plasmodesmata distributions, consistent with biological measurements.
**Plasmodesmata-mediated auxin diffusion accounts for branching patterns in wild-type plants.** With this 3D model, we first aimed to quantitatively assess whether plasmodesmata-mediated auxin diffusion could serve as a plausible inhibitory mechanism regulating branch distribution in *Physcomitrium*. Specifically, we tested whether the model could account for the sizes of the three branch inhibition zones characterizing branching patterns : the apical inhibition zone (AIZ), the inter-branch inhibition zone (IBIZ), and the basal inhibition zone (BIZ).

Previous works have shown that these inhibition zones are influenced by distinct auxin sources (Figure S6)¹³,¹⁵. The BIZ size is likely regulated by auxin derived exclusively from a large source at the gametophore base, while the IBIZ size is thought to depend on the inhibitory strength of the most recent branch. In contrast, the AIZ size appears to be determined by the combined effects of two auxin sources : one at the gametophore apex and the other at the most apical branch. Based on the observed distributions of AIZ and IBIZ size, we estimated the distribution of the AIZ proper (AIZ*) size, representing the inhibitory zone influenced by the apical auxin source only (Figure S6, see STAR Methods). This approach simplified the system's complexity, enabling us to predict the size of each zone by simulating the inhibitory effect of a single, zone-specific auxin source.

To achieve this, we divided the idealized gametophore into independent apical, middle and basal portions, corresponding to the AIZ*, IBIZ or BIZ, respectively (Figure 4A). Plasmodesmata distributions, modeled from longitudinal and transverse histological sections taken at the gametophore apex and at the level of an initiating branch, were applied to the cellular templates representing the apical and middle portions, respectively (Figures 3, 4B and 4C), and *h* was assumed to be constant for the basal portion (Figure 4D). Auxin diffusion through plasmodesmata was then simulated separately for each portion (Figures 4E-4G, see STAR Methods).

The interpretation of auxin concentration across the tissue, which ultimately determines the size of branch inhibition zones, depends on the auxin sensitivity of branch-competent cells. These cells are nested in phyllid axils, with their positions defined by the phyllotactic pattern. Auxin sensitivity was represented by a threshold $τ$ in the model, such that branching was inhibited when the local auxin concentration exceeded a given value of $τ$. Conversely, a branch would initiate when the local auxin concentration in branch-competent cells dropped below that value. Hence, branch inhibition zone sizes were predicted over a 2D space defined by two parameter ratios : $τ/α$, expressing the non-dimensionalized auxin sensitivity threshold as a fraction of the auxin concentration at the source, and $P/δ$, the permeability per plasmodesmata over the degradation rate. Realistic values of $P/δ$ ranged from 300 to 4000 µm³, while $τ/α$ was considered to be above 0.01, accounting for biological variability (Data S1A, Figures S7A-S7C).

To compare the model outputs with biological data, we assumed that $P/δ$ and $τ/α$ varied stochastically according to normal distributions centered on mean values $µ$ with variance $σ^2$. We explored the resulting four-dimensional space defined by ($µ_{P/δ}$, $µ_{τ/α}$, $σ_{P/δ}$, $σ_{τ/α}$), where each parameter set corresponded to a specific distribution of inhibition zone sizes, estimated with 1000 draws. These distributions were then compared to the biological data using the Wasserstein distance, which estimates the cost of transforming one distribution into another (Figures 4H-4J, see STAR Methods). For each gametophore portion, the parameter distribution yielding the best score (i.e., the lowest Wasserstein distance) was shown in the ($P/δ$, $τ/α$) plane (Figures 4K-4M). The comparison of the corresponding inhibition zone size distribution with the data collected from wild-type plants demonstrated that the model successfully captured the observed phenotypes, including both the average and the variability within the interquartile range of the biological data (Figures 4N-4P). Furthermore, we found that a single parameter distribution was also capable of approximately capturing the size range of all three inhibition zones (Figures S7D-S7H). Additionally, extending the model to incorporate growth, along with all auxin sources and inhibition zones at the whole-gametophore level, yielded similar results (Figure S8). This confirmed that modeling separate inhibition zones in fixed tissues does not substantially alter the predictions, while also improving computation time.

Together, these results indicate that a model of branching control based on symplasmic auxin diffusion, calibrated with realistic parameters, consistently predicts branch inhibition zone sizes at the whole gametophore scale. Notably,
the model successfully reproduces the average sizes and much of the variability in branch spacing observed in real plants.

**Changes in the model's plasmodesmal permeability explain branching patterns in plants with altered callose metabolism.** To further assess the plausibility of our model, we tested its capacity to capture the changes in branch patterning caused by alterations in plasmodesmal permeability.

First, to reduce plasmodesmal permeability in *Physcomitrium*, we generated *icals3m* transgenic lines containing an estrogen-inducible version of *cals3m*, a gain-of-function allele of the *Arabidopsis CALS3* gene that enhances callose biosynthesis<sup>21</sup>. We confirmed that *cals3m* induction upon β-estradiol application increased callose deposition at cell-cell junctions and reduced intercellular diffusion in moss, consistent with its expected effects (Figures S9A-F)<sup>34,35</sup>.

To assess the effect of reduced plasmodesmal permeability on gametophore branching, *icals3m* plants were grown for 6-8 weeks, repeatedly soaked in a 1 µM β-estradiol solution, and their branching patterns were quantified as previously described<sup>15</sup>. In control experiments, we noticed that DMSO (0.01% v/v), the solvent used to dissolve β-estradiol, reduced branch production compared to water-treated controls reported in previous studies (see<sup>13,15</sup>), while β-estradiol itself had a mild positive effect on branching in WT plants (Figures S9G-S9K). Despite this, we found that branch density was significantly higher in *icals3m* than in WT following β-estradiol treatment, indicating that reduced symplasmic permeability may promote gametophore branching (Figures S9G-S9K).

Next, we used *icals3m* lines to test the biological relevance of the model's plasmodesmata permeability parameter (*P*). We hypothesized that branching pattern perturbations, specifically a reduction in AIZ and IBIZ size (Figure S9L), could be explained by reduced symplasmic permeability. Consistent with this hypothesis, we found that a reduction in AIZ* and IBIZ size could be captured by the model by decreasing *P* values (Figure 5, Data S1B). However, because the mock treatment reduced branching, apical inhibition zones were larger than usual, making quantitative comparison between biological data and model output within realistic parameter ranges more challenging (see purple values in Figure 5A).

To investigate this regulatory mechanism further, we disrupted endogenous *CALLOSE SYNTHASE (CALS)* genes, which increases symplasmic permeability in other plants, such as *Arabidopsis*<sup>23</sup>. Using our previously published transcriptome dataset<sup>13</sup>, we identified two candidate genes, *CALS4* (*Pp3c10_19330*) and *CALS12* (*Pp3c4_15060*), whose expression is repressed by short-term auxin treatment but induced after gametophore decapitation, suggesting a role in auxin-dependent branching control (Figure S10A). We analyzed *CALS4* and *CALS12* expression patterns using transcriptional fusions between their respective promoters and a β-glucuronidase (GUS) gene introduced into a WT background (Figures 6A-6D, S10B and S10C). We observed overlapping expression domains for *CALS4* and *CALS12* promoters in gametophores, with GUS activity predominantly detected in the basal part of the stem, the main apex, and initiating branch apices. The *CALS12* promoter also exhibited expression in a sub-apical domain of unknown function and in axillary hairs. These observations suggested a probable role for both genes in gametophore branching.

To test this hypothesis, we knocked out *CALS4* and *CALS12* genes using CRISPR/Cas9-mediated genome editing, which resulted in specific deletions of more than 19 and 9 kb, respectively (Figures S10D-S10F). We confirmed that *cals4/12* mutant plants exhibited reduced callose deposition at cell-cell junctions using aniline blue staining (Figures S11A and S11B). To directly assess how these mutations affected intercellular diffusion in gametophores, we developed a novel assay in which Sephadex beads loaded with the diffusible, fluorescent dye HPTS were applied to the cut surface of decapitated shoots. After two hours, diffusion was visualized by confocal microscopy, and the length of the fluorescence signal was quantified. Mutants showed enhanced diffusion compared to the wild type (Figures S11C-S11I). Branching patterns of mutant gametophores were then analyzed after 5-7 weeks of growth. A significant increase of the BIZ size was observed in all mutants, compared to WT controls,
with a more pronounced effect in *cals12* and *cals4/12* backgrounds (Figures S12A-S12H). This result indicates that, while both genes play a role in patterning branch distribution in the proximal region of the gametophore, *CALS12* contributes relatively more significantly.

We hypothesized that mutant phenotypes could be captured by increasing symplasmic permeability in the model. Consistent with this hypothesis, increasing *P* values in the model reproduced the observed increase in BIZ size (Figures 6E-6G, S13; Data S1B).

To investigate further the effect of increasing plasmodesmal diffusivity by reducing callose levels, while limiting possible genetic compensation by other endogenous *CALS* genes, we undertook a gain-of-function approach. To this end, we used a transgenic line that constitutively overexpresses the plasmodesmata-localized callose-degrading protein *PpGHL17-1* (*GHL17-1oe*), leading to increased symplasmic diffusivity<sup>36</sup>. Using qRT-PCR, we confirmed that *GHL17-1* was expressed in gametophores, with higher levels in lateral branches (Figure S11J), validating the suitability of the *GHL17-1oe* line as a tool for modulating cell-to-cell permeability in the tissues under study. We also verified that *GHL17-1oe* plants exhibited reduced callose deposition at cell-cell junctions (Figures S11A and S11B) and enhanced intercellular diffusion, as shown by the HPTS diffusion assay (Figures S11C-S11I). We then analyzed branch distribution in *GHL17-1oe* gametophores and found an overall reduction in branch density, along with a significant increase in IBIZ size, in comparison with WT (Figures S12I-S12L). This branching phenotype could be captured by the model by increasing *P* values (Figures 6H-6J and Data S1B).

Our findings consistently show that increasing callose accumulation promotes branching, while disrupting callose biosynthesis or promoting callose degradation reduces branching. Our model supports the hypothesis that the observed changes in branch distribution in the corresponding transgenic lines may result from restricted or enhanced symplasmic auxin diffusion, respectively.

**3D gametophore geometry and phyllotaxis drive branch spacing robustness.** The above results show that an auxin diffusion-based branching control mechanism can explain much of the size variability in the inhibition zones that define gametophore branching patterns. However, a key unexplained feature of these patterns was the narrow peak in the number of phytomers separating consecutive branches in the IBIZ, centered at 4 (Figure 7A) - a feature that was not seen for the AIZ and BIZ (Figures S6B, S6C and S6E). This peak was indicative of a potentially robust phenotype that could, in principle, either emerge from the diffusion mechanism itself or result from additional, yet unidentified, molecular regulatory processes.

To determine whether this feature could theoretically arise from symplasmic diffusion alone, we used our model to quantify the stability of the predicted outputs under varying input parameters, a key criterion for defining robustness. For this analysis, we defined a robustness index ($R_i$) that represented the proportion of all $\tau/\alpha$ values across the parameter space, for a given $P/\delta$ value, yielding a particular IBIZ size (Figures 7B and S14). We found that the relationship between $R_i$ and IBIZ size was non-monotonic, with maximum $R_i$ values consistently observed at 3-4 phytomers, regardless of the $P/\delta$ value, reflecting a robust model outcome. This theoretical result supports the idea that the observed peak in IBIZ size distribution is likely an emergent property of the diffusion mechanism alone, and that an as-yet unidentified feature constrains the model output under variable input parameters.

To investigate further the source of IBIZ robustness, we reasoned that a major factor differentiating the IBIZ from the two other inhibition zones was the distance separating the corresponding auxin sources from the successive phyllid axils. Considering that symplasmic diffusivity is more limited by the number of cells than by physical distance, we defined this distance as the length - measured in cell number - of the shortest path between an auxin source and a phyllid axil within the cellular grid of the 3D gametophore. This distance increases linearly for apical and basal auxin sources (Figures S15A and S15B), whereas it oscillates for a peripheral auxin source (Figure 7C and 7D). This difference is explained by the radial component of the distance, which varies between 0 and the whole stem diameter in the case of a peripheral auxin source, depending on phyllotaxis. In particular, this results
in a local maximum at the fourth phyllid from the auxin source (black arrow in Figure 7C). These observations suggested that spiral phyllotaxis imposes a constraint, making some phyllid axils more susceptible to inhibition than others based on their position relative to the closest peripheral auxin source. Specifically, the fourth phyllid axil from the source is less prone to auxin inhibition than the third or the fifth. We therefore hypothesized that this constraint underlies the observed robustness of the IBIZ at 4 phytomers. As it depends on the relative contributions of radial and axial distances between phyllid axils and the peripheral auxin source, we reasoned that changes in stem diameter or phytomer length could modulate IBIZ robustness (Figures S15C and S15D).

To test this, we simulated auxin diffusion in the middle section of idealized 3D gametophores with spiral phyllotaxis, varying stem diameter and phytomer length, and computed the robustness index $R_i$ for IBIZ sizes of 2-10 phytomers (Figures 7E-7G). Simulation outputs were compared with a control situation where divergence angles were null, *ie.* with phyllids aligned vertically above one another (Figures 7H-7J and S15E-S15G). We observed that reducing phytomer length (Figure 7F) or increasing stem diameter (Figure 7G) in gametophores with spiral phyllotaxis progressively increased robustness for an IBIZ size of 4 phytomers. However, this effect was not observed in gametophores with aligned phyllids (Figures 7I and 7J), supporting our hypothesis.

This phenomenon suggested additional implications for branching regulation. With radial diffusion emerging as a critical factor, we hypothesized that the relative plasmodesmata densities on longitudinal *versus* transverse cell walls would significantly contribute to IBIZ robustness. To verify this, we simulated a scenario with uniform plasmodesmata densities on longitudinal and transverse cell walls, effectively representing the absence of plasmodesmata dilution on longitudinal walls. This change increased radial symplasmic diffusivity in the gametophore stem, which abolished IBIZ size robustness at 4 phytomers (Figures S15H-S15L).

A key emergent property of the symplasmic diffusion mechanism predicted by our model was the increase in IBIZ size robustness with stem diameter (Figure 7G), which we sought to investigate experimentally. Although mosses lack dedicated stem thickening tissues like the cambium, we anticipated that gametophore diameter would increase with age as a result of cell division occurring throughout the stem and, consequently, that stem thickening would enhance IBIZ size robustness during development. To verify this, we first plotted stem diameter against length (used as a proxy for age) in WT plants and confirmed that gametophore stems gradually increase in thickness over time (Figure 7K). Next, using our model and previously optimized parameter distributions capturing WT branching patterns (Figure 4L), we predicted IBIZ size distributions for various stem diameters (Figures 7L-7O). Simulations indicated that IBIZ size would be more variable with smaller stem diameters but would become increasingly stable, converging towards 4 phytomers as stem diameter increases (Figure 7P). We tested this prediction by examining the spatio-temporal variability of IBIZ size in real plants, using the relative position of branches as a proxy for their age and developmental sequence. This approach relies on the assumption that older branches, near the gametophore base, develop when the stem is narrower, whereas younger branches, located closer to the gametophore tip, emerge when the stem is comparatively wider¹⁵. IBIZ size distribution was therefore plotted according to the rank along the gametophore axis, considering only gametophores with three or more branches (Figures 7Q and 7R). As anticipated, IBIZ size was more variable in rank 1, corresponding to the earliest-formed branches, and became more consistent, converging at 4 phytomers with increasing ranks. These measurements validated the prediction that stem thickening enhances IBIZ size robustness during gametophore development.

Together, these results demonstrate that robust branch spacing can emerge from symplasmic auxin diffusion through its interplay with 3D gametophore geometry and spiral phyllotaxis. They also highlight the critical role of limiting radial symplasmic diffusivity in maintaining IBIZ robustness.
Discussion

This work aimed to explore the mechanistic basis of branch patterning, a key determinant of plant architecture, in the moss *Physcomitrium patens*. We have taken a 3D computational modeling approach to establish a comprehensive framework integrating key aspects of auxin physiology, plasmodesmata distribution and function, and gametophore geometry, effectively bridging molecular, cellular and organ scales, and based on realistic biological parameters. Our results demonstrate that branching control by auxin-mediated inhibition, relying solely on symplasmic diffusion for intercellular auxin movement, can accurately account for observed branching patterns at the whole gametophore level.

Mapping plasmodesmata in wild-type gametophores has revealed graded distribution patterns at the tissue level, with the highest densities in meristematic regions. We observed that plasmodesmal connectivity decreases with tissue differentiation and is lower on longitudinal walls than transverse walls in mature tissues, likely due to cell wall expansion, which suggests a greater diffusion capacity along the apical-basal axis than in the radial orientation. Consistent with this, *Physcomitrium* phyllids principally contain primary plasmodesmata, and secondary plasmodesmata do not form at a sufficient rate to compensate for wall elongation³⁷. These findings also reveal shared key features with the morphologically distinct meristems of angiosperms and ferns³⁸, where plasmodesmata densities are generally higher than in surrounding tissues, cell differentiation is associated with reduced plasmodesmata number and cell-to-cell diffusion, and molecular transport primarily occurs in the longitudinal orientation²⁶,³⁹⁻⁴⁶.

Besides plasmodesmata densities, previous work in vascular plants has shown that plasmodesmal clustering can reduce effective symplasmic permeability⁴⁷. Although the extent of clustering in *Physcomitrium* gametophores remains unknown, our electron microscopy analyses did not reveal evidence of such clustering. Moreover, plasmodesmal clustering is typically associated with the formation of secondary plasmodesmata⁴⁸. Since secondary plasmodesmata are rare in *Physcomitrium*³⁷, it is reasonable to assume that clustering is limited in gametophore tissues and unlikely to significantly affect symplasmic permeability. The structural and functional properties of individual plasmodesmata can also affect the diffusion capacity of the symplasmic network. These properties can vary across organs and are dynamically regulated over time⁴⁹. Detailed data are now available at the scale of a few cells⁴⁷,⁵⁰, suggesting that incorporating this variability might enhance the model’s predictive accuracy. However, such data remain inaccessible at larger scales, and collecting them is beyond the scope of this study. Moreover, adding more data into a model and increasing its complexity does not necessarily improve its accuracy⁵¹. It is important to stress that the benefit of our approach lies in integrating key averaged parameters influencing symplasmic diffusivity within a unified framework. This enabled us to effectively connect different scales and make predictions at the level of an entire shoot, which would be extremely challenging, if not impossible, to achieve otherwise.

Studies in Arabidopsis have established that auxin diffuses through the symplasmic pathway, regulated by callose-mediated plasmodesmal gating¹⁸. For instance, the *Arabidopsis CALS10/GSL8* gene is directly activated by *ARF7*-dependent auxin signaling, creating a feedback loop where auxin modulates its own diffusion²³. We have found that short-term auxin treatment represses *CALS4* and *CALS12* gene transcription in *Physcomitrium*, suggesting that a conserved auxin-callose interplay has been remodeled during land plant evolution. Nevertheless, it remains unclear whether the effect of auxin persists over time, and how the spatiotemporal dynamics of this feedback influence diffusion in the long term. Regardless, CRISPR/Cas9-mediated deletion of these genes resulted in reduced gametophore branching, a phenotype that our model captured by increasing plasmodesmal permeability and auxin diffusion. The phenotype was restricted to the gametophore base, indicating possible genetic redundancy. To bypass this, we took gain-of-function approaches, promoting callose degradation through ectopic *GHL17-1* expression and enhancing callose accumulation via conditional *cals3m* activation. Both strategies resulted in more pronounced branching perturbations than those observed in *cals4/12* mutants, consistent with our working hypothesis.

The absence of major developmental defects and the differing branching phenotypes across inhibition zones in transgenic lines suggest that callose homeostasis is tightly regulated, with compensatory mechanisms likely pre-
venting detrimental effects. Notably, these compensatory mechanisms may act heterogeneously within individual gametophores. Indeed, we found no correlation between apical and inter-branch inhibition zone sizes within individual gametophores, in both transgenic lines and wild-type plants (Figure S6D). If variation in zone size primarily reflected inter-individual differences in symplasmic permeability (or related parameters), one would expect these zones to correlate. The lack of correlation instead indicates intra-individual variability, implying that regulatory mechanisms influencing branch inhibition are spatially heterogeneous within a single gametophore.

Nevertheless, it is also possible that, in specific regions of the gametophore, regulatory mechanisms other than callose-dependent pathways may predominate<sup>49,52</sup>. For instance, an intercellular flow control mechanism involving Multiple C2 domains and Transmembrane region Proteins (MCTPs), which tether the desmotubular endoplasmic reticulum to the plasma membrane in plasmodesmata<sup>53</sup>, has recently been shown to override callose effects in *Arabidopsis*<sup>54</sup>. While MCTPs are evolutionary conserved plasmodesmal proteins<sup>36,55</sup>, their role in *Physcomitrium* remains unclear, making them prime candidates for investigating alternative regulatory mechanisms of symplasmic diffusion.

Beyond auxin, other phytohormones, proteins and RNAs may also traffic through plasmodesmata<sup>49,56,57</sup>. While our model convincingly demonstrates that simulating auxin diffusion alone is sufficient to capture wild-type branching patterns and phenotypic changes in transgenic lines with altered callose biosynthesis, it cannot be ruled out that the movement of other molecules is not affected in those lines. For example, unknown plasmodesmata-mobile branch-promoting signals might counterbalance auxin's effect, possibly explaining phenotypic disparities in transgenic gametophores. A key challenge therefore remains to identify alternate mobile factors that can exert control over branch outgrowth, if any exist, potentially via targeted genetic screens similar to those investigating KNOT-TED1 transcription factor movement<sup>58</sup>.

Molecular diffusion is typically viewed as a slow process, sensitive to noise within the cellular environment, generating short-range gradients over only a few cells<sup>59</sup>, which makes it seemingly incompatible with the formation of robust developmental patterns. Our modeling experiments indicate that symplasmic diffusion can yield molecular gradients reaching distances equivalent to observed branch inhibition zone sizes within biologically realistic parameter value ranges. Regarding the IBIZ size, we have found that its variability remained low around a dominant value of four phytomers and decreased during development, eventually becoming nearly negligible, which initially seemed inconsistent with a diffusion-based lateral inhibition mechanism. 3D simulations in an idealized gametophore have resolved this paradox, demonstrating that robust branching patterns can emerge from the interplay between symplasmic diffusion, organ and tissue geometry, and phyllotaxis. Notably, the model predicts that stem thickening may reduce IBIZ size variability over time. It also highlights the important role of phyllotaxis as a pre-pattern that constrains potential branch initiation sites along the gametophore axis.

We conclude that symplasmic auxin diffusion may play a more significant role in regulating shoot branching in *Physcomitrium* compared to vascular plants like *Arabidopsis* or *Selaginella*, where PIN-mediated polar auxin transport is predominant<sup>10,16</sup>. The reasons why PIN proteins do not mediate long-range polar auxin transport in the stem of moss gametophores remain unclear. Nonetheless, it is important to note that PINs are important for meristem function in *Physcomitrium*, likely establishing auxin gradients around the apical cell<sup>14,60</sup>. Given that these diffusive and active transport mechanisms have likely co-existed since the emergence of land plants<sup>61</sup>, we speculate that the observed differences between species may be more related to their size than their phylogenetic position. For instance, *Physcomitrium* gametophores are typically about half a centimeter tall, whereas *Arabidopsis* inflorescence shoots can reach up to 50 cm, representing a difference of two orders of magnitude. This supports the idea that distinct auxin transport capacities are required in each species, with long-range transport being crucial in the latter situation. Interestingly, other moss species can have gametophores that exceed dozen centimeters in size and display complex branching architectures, unlike *Physcomitrium*<sup>62</sup>. A major challenge moving forward will be to decipher how the various regulatory mechanisms of auxin movement have evolved and contributed to the
diversification of plant architecture.

## Resource Availability

### Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Yoan Coudert (yoan.coudert@cnrs.fr).

## Materials Availability

Plasmids and plant materials generated in this study are all available from the Lead Contact upon request.

## Data and Code Availability

The code developed for this study is available at https://gitlab.inria.fr/mosaic/publications/abitbol-spangaro et_al. Numerical results and experimental data corresponding to auxin degradation assays (Fig. S1), plasmodesmata maps (Fig. S3), Dendra2 diffusion assays (Fig. S9), and branching patterns of wild-type gametophores (Fig. S6) and *cals4-12* mutant gametophores (Fig. S12) can be found on the same public repository. Further requests for resources should be directed to the Lead Contact.

## Acknowledgements

Y.C. acknowledges Jill Harrison for her support at the onset of this project, as well as Mitsuyasu Hasebe for the pPGX plasmid, Tomomichi Fujita for the Dendra2 plasmid, Ykä Helariutta for the *cals3m* gene, Jan Petrášek for discussions about auxin degradation, Rüdiger Simon for the *GHL17-1oe* line and Jonathan Rudolph for his assistance with the qPCR experiment. J.A.-S. acknowledges Guillaume Cerutti for assistance with the cellcomplex Python library. J.d.K. acknowledges Tjitske Riksen for assistance with intercellular diffusion measurements. Research in Y.C.'s team was financially supported by the CNRS (ATIP-Avenir grant, Ph.D. fellowship to J.A.-S.), the ANR (grant ANR-21-CE20-0024), and the Museum National d'Histoire Naturelle (ATM grant). Plasmodesmata imaging was performed at the Imagerie-Gif core facility, supported by the Agence Nationale de la Recherche (ANR-10-INBS-04/FranceBioImaging; ANR-11-IDEX-0003-02/ Saclay Plant Sciences). V.C. is a Research Associate of the Fonds de la Recherche Scientifique - FNRS and was funded under the grant number 1208619F at the onset of the project.

## Author contributions

Conceptualization: J. A.-S., C. Go. and Y.C.; Formal analysis and software: J. A.-S.; Funding acquisition: Y.C.; Investigation and resources: J. A.-S., G.C., A.M., S.H., C.B., C.Gi., V.S., P.I.D., R.S., V.C., J.d.K. and Y.C.; Project administration: Y.C.; Supervision: C. Go. and Y.C.; Visualization: J. A.-S. and Y.C.; Writing: J. A.-S. and Y.C., with input from C.B., V.S., P.I.D., R.S., V.C., J.d.K. and C. Go.

## Declaration of interests

The authors declare no competing interests.
Legends for Main Figures

**Figure 1.** Modeling symplasmic auxin diffusion in an idealized gametophore. (A-B) Wild-type gametophore before (A) and after (B) removing the phyllids of the primary stem. Asterisks indicate lateral branches. Scale bar = 1 mm. (C) 3D representation of the idealized gametophore tissue created by extruding cross-sections generated through Voronoi tessellations. (D) Schematic representation of the auxin diffusion model parameters, as defined in Eq.1. Auxin concentration *a* is computed in each cell. Auxin biosynthesis is modeled by setting a constant concentration *α* in a set of cells *S* (shown in red). For any cell *i* in the tissue (shown in light blue), the temporal variation in auxin concentration is determined by the sum of the diffusive fluxes with each neighboring cell *j*. The diffusive flux between cells *i* and *j* is proportional to the cell wall permeability, which depends on the plasmodesmata density on the corresponding cell wall ρ(*i*, *j*), the plasmodesmata open cross-section *ω*, and the cell wall thickness *l*. Abbreviations: DT, desmotubule; PD, plasmodesma. See also Figure S1 and Data S1A.

**Figure 2.** Plasmodesmata exhibit a graded distribution in gametophores, with maximal densities in meristematic regions. (A) Schematic diagram of a moss gametophore with a lateral branch. Lines indicate approximate locations of sectioning. (B-D) Representative transverse (B, C) and longitudinal (D) sections of moss gametophores, with plasmodesmata densities mapped onto corresponding walls. Arrowheads mark apical cells of a lateral branch (C) and a primary gametophore axis (D). (E, G) Box plots showing that plasmodesmata densities are overall higher at the main gametophore apex than in other tissues (E), and on transverse walls than longitudinal walls (G). (F, H) Dot plots showing the decreasing gradient of plasmodesmata densities as a function of the distance from secondary and primary apical cells in sections (C) and (D), respectively. Scale bars = 10 µm (B) and 50 µm (C, D). See also Figures S2 and S3.

**Figure 3.** Modeling observed plasmodesmata densities in idealized gametophore tissues. (A-B) Projection of estimated plasmodesmata densities onto idealized transverse (A) and longitudinal (B) tissues. Plasmodesmata densities depend on the distance to an auxin source, defined as the epidermal cell marked by an arrowhead in (A) and the entire apical cell layer indicated by brackets in (B). (C-E) Purple curves represent simulated plasmodesmata densities in idealized tissue sections corresponding to the transverse (C, longitudinal walls) and longitudinal (D, transverse walls; E, longitudinal walls) orientations. These curves are overlaid with dots representing plasmodesmata densities quantified in corresponding real tissues (Figures 2C and 2D). Scale bars = 20 µm. See also Figures S4 and S5.

**Figure 4.** A model based on symplasmic auxin diffusion captures branch inhibition zone sizes in Physcomitrium gametophores. (A) Apical, middle and basal portions of an idealized 3D gametophore used to simulate the apical inhibition zone proper (AIZ*), inter-branch inhibition zone (IBIZ) and basal inhibition zone (BIZ). Phytomers are depicted in alternating shades of grey, with branch-competent cells in leaf axils shown in light green. Auxin sources are marked in red. Red frames indicate the gametophore's regions shown in panels B-G. (B-G) Selected regions of the gametophore's apical (B, E), middle (C, F) and basal (D, G) portions, displaying modeled plasmodesmata densities projected onto transverse and longitudinal cell walls, with auxin sources in red (B-D), and simulated auxin concentrations at equilibrium in individual cells (E-G). (H-P) Comparison of model simulations with biological data for the AIZ* (H, K, N), IBIZ (I, L, O) and BIZ (J, M, P). (H-J) Heat maps of the Wasserstein distance, representing the similarity between the predicted and observed phenotype distributions in wild-type (WT) plants. Heat maps show minimum distances projected on the plane ($\mu_{P/\delta}, \mu_{\tau/\alpha}$) from the four-dimensional space ($\mu_{P/\delta}, \mu_{\tau/\alpha}, \sigma_{P/\delta}, \sigma_{\tau/\alpha}$). White stars indicate the parameter sets best matching the biological data. (K-M) Scatter plots of 1000 samples (white dots) drawn from normal distributions corresponding to white stars in panels (H-J) on the parameter space ($P/\delta, \tau/\alpha$). Colors indicate the predicted inhibition zone sizes for
given parameter values. (N-P) Boxplots showing that inhibition zone size distributions predicted by the model, as depicted in panels (K-M), align well with the observed measurements in WT plants. See also Figures S6-S8 and Data S1B.

**Figure 5. Decreasing plasmodesmal permeability in the model captures branching pattern perturbations in *icals3m* gametophores after induction of callose deposition.** Comparison of the model output with biological data for the apical inhibition zone proper (AIZ\*) (A,C,E) and inter-branch inhibition zone (IBIZ) (B,D,F) in *icals3m* gametophores following mock or β-estradiol treatment. (A-B) Heatmaps showing the top 2.5% of parameter values that minimize the distance between model predictions and observed sizes of branch inhibition zones in *icals3m* gametophores after mock (purple) or β-estradiol (fawn) treatment. (C-D) Scatter plots of 1000 samples drawn from normal distributions yielding the highest similarity between simulated and biological data in the parameter space ($\mu_{P/\delta}, \mu_{\tau/\alpha}$), with distributions corresponding to mock-treated plants in white and β-estradiol-treated plants in brown. Predicted inhibition zone sizes corresponding to specific parameter values are displayed on a color range from dark blue to light grey. (E-F) Boxplots showing that, while high values of AIZ\* size in mock-treated gametophores can not be reproduced by the model within realistic parameter values, reducing the parameter *P* corresponding to plasmodesmal permeability captures the decrease in AIZ\* (E) and IBIZ (F) sizes observed between mock-treated and β-estradiol-treated plants. See also Figure S9 and Data S1B.

**Figure 6. Increasing plasmodesmal permeability in the model captures branching pattern perturbations in *cals4/12* mutants and *PpGHL17-1* overexpressors. (A-D) GUS staining of *CALS4* (A, B) and *CALS12* (C, D) transcriptional reporter lines showing that *CALS4* and *CALS12* promoter activity in overlapping domains, including the main apex (black arrowheads) and the middle-basal region (brackets) of gametophores (A, C). GUS activity was also detected in emerging lateral branches, including rhizoids (black arrows) for both genes, and in axillary hairs (white arrowheads) for *CALS4* (B, D). Dashed lines mark the boundary between the stem and the detached phyllid (B, D). (E-G) Comparisons between model simulations and biological measurements of basal inhibition zone (BIZ) size in *cals4/12* double mutants and wild-type (WT) controls. (E) Heatmaps showing the top 2.5% of parameter values that minimize the distance between model predictions and biological data for *cals4/12* mutants (fawn) and WT (purple). (F) Scatter plot of 1000 samples drawn from normal distributions that maximize the similarity between simulated and biological data in the parameter space ($\mu_{P/\delta}, \mu_{\tau/\alpha}$), with WT data in white and *cals4/12* mutant data in brown. Colors indicate the predicted inhibition zone sizes for corresponding parameter values. (G) Boxplots showing that the distribution of BIZ sizes predicted by the model, as represented in panel (F), align with the biological data from WT and *cals4/12* gametophores. (H-J) Comparison between model simulations and biological measurements of inter-branch inhibition zone (IBIZ) size in *PpGHL17-1* overexpressors (*GHL17-1oe*) and WT controls, using similar data representation as in (E-G). The distribution of IBIZ sizes predicted by the model matches the biological data from WT and *GHL17-1oe* gametophores. Scale bars = 1 mm (A, C), 100 $\mu$m (B, D). See also Figures S10-S13, and Data S1B and S1C.

**Figure 7. 3D gametophore tissue geometry and phyllotaxis underpin branch spacing robustness.** (A) Histogram representing the frequency of observed IBIZ sizes in wild-type gametophores, with a peak at 4 phytomers. (B) Plot of the robustness index $R_i$ for IBIZ sizes predicted by the model, showing maximal $R_i$ values at 4 phytomers (black arrow) for three different $P/\delta$ values. (C) Plot showing that, in a gametophore with spiral phyllotaxis (137.5°), the distance in cell numbers of consecutive phyllids from an axillary auxin source (purple line) is determined by a linear axial component (dotted blue line) and an oscillating radial component (dashed red line), resulting in a local maximum at the 4^th^ phyllid (black arrow). (D) 3D representation of the distance of the 4^th^ phyllid from an axillary auxin source in an idealized 3D gametophore. Only cell layers are represented, with branch-competent cells in leaf axils shown in green, and the auxin source in red. Axial, radial and total distances
follow the same color coding as in (C). (E) Phyllid arrangement in a gametophore with spiral phyllotaxis. Numbers indicate the phyllid index (green) relative to the auxin source (red). (F-G) Heatmaps showing that IBIZ sizes corresponding to maximum $R_i$ values (white dots) converge toward 4 phytomers with decreasing phytomer length (F) and increasing stem diameter (G). $R_i$ values are displayed using a green-to-yellow color scale. Stem diameter is fixed at 14 cells in (F), and phytomer length is fixed at 3 cell layers in (G). (H) Phyllid arrangement in a gametophore with null divergent angles, represented similarly to (E). (I-J) Heatmaps showing that, for aligned phyllids, $R_i$ decreases monotonically with IBIZ size regardless of phytomer length (I) or stem diameter (J). Data representation is consistent with (F-G). (K) Boxplots showing that stem diameter in wild-type gametophores gradually increases during development, ranging from 150 to 380µm. Data are grouped by gametophore height ($n = 12, 11, 18, 14, 14$) and position along the gametophore axis, with three measurements taken per apical, middle, and basal region per gametophore. Dotted lines correspond to simulated stem diameters shown in (L-P). (L-M) Cross-sections of idealized gametophores with diameters, denoted $\emptyset$, of 12 (L) and 16 (M) cells, corresponding to 240 and 320 µm, respectively. Scale bars = 20 µm. (N-O) Scatter plots of 1000 samples (white dots) from a normal distribution optimized to reproduce wild-type branching patterns (Figure 5L) in the parameter space ($P/\delta, \tau/\alpha$). Colors indicate predicted inhibition zone sizes for idealized gametophores with diameters of 12 (N) and 16 (O) cells. (P) Boxplots showing that IBIZ size variability predicted by the model for stem diameters of 12, 14 and 16 cells decreases and converges towards 4 phytomers with increasing stem diameter. (Q) The IBIZ rank $n$ is defined as the interval between the $n^{th}$ and $n + 1^{th}$ branches from the base. (R) Boxplot showing that IBIZ size variability decreases with IBIZ rank in wild-type gametophores. Data match the model predictions (P), suggesting that stem thickening during gametophore development enhances IBIZ robustness at 4 phytomers. Phytomer and phyllid numbering is consistent across all panels. See also Figures S14 and S15.
The Gransden strain of *Physcomitrium patens* was used as the wild-type (WT) background for generating all the transgenic lines. Moss colonies were initiated from protonemal cultures and cultivated on BCDAT medium (250 mg/L MgSO₄·7H₂O, 250 mg/L KH₂PO₄ (pH 6.5), 1010 mg/L KNO₃, 12.5 mg/L, FeSO₄·7H₂O, 920 mg/L di-ammonium tartrate (C₄H₁₂N₂O₆), and 0.001% Trace Element Solution (0.614 mg/L H₃BO₃, 0.055 mg/L AlK(SO₄)₂·12H₂O, 0.055 mg/L CuSO₄·5H₂O, 0.028 mg/L KBr, 0.028 mg/L LiCl, 0.389 mg/L MnCl₂·4H₂O, 0.055 mg/L CoCl₂·6H₂O, 0.055 mg/L ZnSO₄·7H₂O, 0.028 mg/L KI and 0.028 mg/L SnCl₂·2H₂O)), supplemented with 7 g/L agar and CaCl₂ added to a 1 mM concentration after autoclaving. Standard growth conditions were 23°C under a 16h light/8h dark cycle, at 50-150 µmol.m⁻².s⁻¹ in MLR-352 growth cabinets.

## Method Details

### Experimental protocols

### Molecular cloning and transgenic line production

**Gene identifiers** Genomic sequences of *PpCALS4* (Pp3c10_19330) and *PpCALS12* (Pp3c4_15060) genes were obtained from the version 3.3 of *P. patens* genome annotation available on Phytozome (https://phytozome-next.jgi.doe.gov)⁶³. The *icals3m* gene is a mutated form of *AtCALS3* (At5g13000), whose genomic sequence is available on TAIR (https://www.arabidopsis.org)²¹. *PpGHL17-1* gene identifier is Pp3c17_13760.

**Inducible *cals3m* (*icals3m*) overexpression** To construct the pPGX6-cals3m expression vector, a 9206 bp DNA fragment corresponding to the *cals3m* genomic sequence, previously cloned into the pDONR-221 plasmid (a gift from Prof. Ykä Helariutta's laboratory), was transferred into the pPGX6 vector (a gift from Prof. Mitsuyasu Hasebe's laboratory) using the Gateway LR Clonase II enzyme. The resulting plasmid was named pPGX6-icals3m and was subsequently transformed into wild-type *P. patens*.

**Photoconvertible Dendra2 constitutive expression** The previously selected *icals3m* line was genetically transformed with the pT10G-Dendra2 plasmid (a gift from Prof. Tomomichi Fujita), which drives the expression of the photoconvertible Dendra2 fluorescent protein under control of the *PpEF1α* promoter.

**PpCALS4 and PpCALS12 promoter reporters** To construct *CALS* promoter reporter vectors, a modified version of the pPIG1-NGGII vector⁶⁴ was first generated. The pPIG1-NGGII vector was digested with ClaI and BlpI to remove the Blasticidin resistance cassette, which was replaced with a Hygromycin resistance cassette previously assembled by PCR, yielding the vector pDEST-MOSSY1. In parallel, an R1/R2 Gateway cassette was amplified by PCR from pT2N2x35SGate⁶⁵, subcloned into pGEMT-easy, and cloned in the correct orientation into the XbaI site of the pDEST-MOSSY1 vector, yielding the pDEST-MOSSY2 vector. A 2743 bp genomic fragment upstream of the *PpCALS4* coding sequence, referred to as the *PpCALS4* promoter, and a 2528 bp genomic fragment upstream of the *PpCALS12* coding sequence, referred to as the *PpCALS12* promoter, were amplified by PCR with the primer pairs AM43F/AM8R and AM10F/AM10R, respectively. Both fragments were then subcloned into the pDONR-221 vector using the Gateway BP Clonase II enzyme and transferred to the pDEST-MOSSY2 vector using
the Gateway LR Clonase II enzyme, yielding the plasmids pDEST-MOSSY2-CALS4pro and pDEST-MOSSY2-CALS12pro, respectively. Both plasmids were linearized with the PmeI enzyme prior to genetic transformation into *P. patens*.

**CRISPR/Cas9-mediated cals deletion mutants** The CRISPR-Cas9 vector system developed by Magdalena Bezanilla's laboratory<sup>66</sup> was used to generate *PpCALS4* and *PpCALS12* knock-out constructs. A series of protospacers was designed using CRISPOR software (http://crispor.tefor.net). These protospacers corresponded to short guide RNAs (sgRNAs) targeting exon 1 (position 130), exon 4 (position 77), exon 33 (position 73) and exon 49 (position 53) of the *PpCALS4* genomic sequence, and exon 2 (position 52), exon 6 (position 75), exon 23 (position 212) and exon 26 (position 124) of the *PpCALS12* genomic sequence. Corresponding oligonucleotide pairs (sgR-PpCALS4-EX01-130-F and sgR-PpCALS4-EX01-130-R, sgR-PpCALS4-EX04-77-F and sgR-PpCALS4-EX04-77-R, sgR-PpCALS4-EX33-73-F and sgR-PpCALS4-EX33-73-R, and sgR-PpCALS4-EX49-53-F and sgR-PpCALS4-EX49-53-R for *PpCALS4*, and sgR-PpCALS12-EX02-52-F and sgR-PpCALS12-EX02-52-R, sgR-PpCALS12-EX06-75-F and sgR-PpCALS12-EX06-75-R, sgR-PpCALS12-EX23-212-F and sgR-PpCALS12-EX23-212-R, and sgR-PpCALS12-EX26-124-F and sgR-PpCALS12-EX26-124-R for *PpCALS12*) were annealed and inserted downstream of a U6 promoter, into the *Bsa*I site of the pENTR-PpU6P-sgRNA-L1R5, pENTR-PpU6P-sgRNA-L5L4, pENTR-PpU6P-sgRNA-R4R3 or pENTR-PpU6P-sgRNA-L3L2 entry vectors (a gift from Magdalena Bezanilla). The different fragments were then transferred using the Gateway LR Clonase II enzyme to either the pMH-Cas9-gate or the pMK-Cas9-gate destination vectors, which contain the Cas9 enzyme coding sequence and confer hygromycin or geneticin resistance, respectively. The resulting plasmids were named pMH-Cas9-PpCALS4 and pMK-Cas9-PpCALS12.

**P. patens genetic transformation** Plasmid DNA was purified using the NucleoBond Xtra Maxi kit. For stable transgene integration, DNA was linearized with appropriate restriction enzymes before transformation into moss protoplasts. Five-day-old protonemal tissue, sub-cultivated at least twice by blending and grown on BCDAT plates overlaid with autoclaved cellophane disks, was used for protoplast preparation. The protonemal tissue was digested in a solution containing 1% driselase and 8% mannitol solution for 60-75 minutes at room temperature in the dark. The resulting protoplast suspension was filtered through a 100 µm nylon mesh, and cells were pelleted by centrifugation at 90 g for 4 minutes, then washed twice in an 8% mannitol solution. Density of the protoplast solution was counted using a Neubauer hemocytometer. Protoplasts were pelleted by centrifugation and re-suspended at a density of 3.2 million cells/ml in an 8% mannitol solution. A 150 µl aliquot of the protoplast suspension was mixed with 150 µl 2xMaMg solution (61 g/L MgCl₂, 0.92 g/L di-ammonium tartrate (C₄H₁₂N₂O₆), mannitol 80 g/L, MES 2 g/L, pH 5.6) and 300 µl PEG solution (23.6 g/L Ca(NO₃)·4H₂O, 4.76 g/L HEPES, 72.8 g/L mannitol, 400 g/L PEG6000, pH 7.5) in a 15 mL tube. The mixture was heat shocked at 45°C for 5 minutes, cooled down at 20°C for 10 minutes and progressively diluted over one hour with 8% mannitol to reach a final volume of 15 mL. Protoplasts were then centrifugated at 90 g for 4 minutes, resuspended in PRM-L solution (BCDAT supplemented with 8% mannitol and 10 mM CaCl₂), and incubated overnight at 23°C in the dark. The next day, protoplasts were centrifugated at 90 g for 4 minutes and resuspended in PRM-T (BCDAT supplemented with 8% mannitol, 10 mM CaCl₂ and 0.4% agar). Protoplasts were plated onto PRM-B solid medium (BCDAT supplemented with 8% mannitol, 10 mM CaCl₂ and 0.8% agar) overlaid with autoclaved cellophane disks, and incubated under standard growth conditions described above. After 5 days, cellophane disks were transferred to BCDAT medium supplemented with antibiotic for selection of transformed protoplasts (50 mg/ml geneticin, 50 mg/ml hygromycin B, or 100 mg/ml zeocin). After 7 days, cellophane disks were transferred to BCDAT medium without antibiotics. A second round of antibiotic selection was performed 14 days later for transformations involving stable insertion of the antibiotic resistance cassette (i.e., *icals3m*, *icals3m Dendra2*, *PpCALS4* and *PpCALS12* promoter reporters). For
CRISPR/Cas9 approaches involving transient selection, only one round of antibiotic selection was performed. All resistant lines were subjected to genotyping by PCR, followed by RT-PCR and/or sequencing analyses as described elsewhere.

## Genotyping and selection of transgenic lines

**PCR** Genotyping was carried out by PCR using chloroform-isoamyl alcohol-extracted gDNA as a template. GoTaq G2 DNA polymerase was used for PCR amplification. *PpUBI* (Pp3c12.5760V3.1) was amplified in all samples using the primer pair YC-PpUBI-152F and YC-PpUBI-152R to confirm DNA quality. Standard thermocycling conditions were: 94°C for 2 min; 94°C for 30 sec, 55-58°C for 30 sec, 72°C for 60 sec/kb (30-35 cycles); 72°C for 5 min. Primer melting temperatures were between 58 and 60°C.

For the putative *cals4*, *cals12* and *cals4/cals12* mutants, clones selected for their antibiotic resistance were genotyped using the primer pairs PD40_F1 and PD40_R1, PD40_F2 and PD40_R2, PD40_F3 and PD40_R3, and PD40_F4 and PD40_R4 for *PpCALS4*; and PD40_F9 and PD40_R9, PD40_F10 and PD40_R10, PD40_F11 and PD40_R11, PD40_F12 and PD40_R12, and PD40_F13 and PD40_R13 for *PpCALS12*. Transgenic lines harboring full deletions of *PpCALS4* and/or *PpCALS12* genomic fragments were sent for Sanger sequencing to confirm the nature of the deletion and predict the resulting protein sequence. Lines with the most drastic mutations were retained for phenotypic analyses.

Predicted putative off-target genomic sites were checked by PCR using the primer pairs PD58-chr21-F and PD58-chr21-R, PD58-chr22-F and PD58-chr22-R, PD58-chr4-F and PD58-chr4-R, PD46_CALS12-Chr27_F and PD46_CALS12-Chr27_R, PD46_CALS12-Chr13_F and PD46_CALS12-Chr13_R, PD46_CALS12-Chr12_F and PD46_CALS12-Chr12_R, PD46_CALS12-Chr08_F and PD46_CALS12-Chr08_R. All putative off-targets had a wild-type sequence, indicating a lack of non-specific mutations at those sites.

Promoter reporter lines were screened by PCR for correct transgene integration into the PIG1 locus using the primer pairs YCL-152F-AM113F and YCL-153R-AM114R, and YCL-153F-AM114F and YCL152R-AM113R. Lines showing correct transgene integration were retained for GUS assays.

For the inducible *cals3m* lines, antibiotic resistant clones were screened by RT-PCR for their response to 1 µM β-estradiol. Clones showing no or reduced *cals3m* expression in response to the solvent control DMSO and a strong induction following β-estradiol treatment were retained for further experiments.

For the *icals3m* Dendra2 lines, following selection based on zeocin resistance, the lines were first screened by PCR with the primer pair YCL-58F and YCL-58R to check for the integration of Dendra2 coding sequence. Subsequently, their basal Dendra2 expression levels were measured, and the lines displaying the brightest fluorescence levels were retained for further analyses. A wild-type line carrying the pT10G-Dendra2 vector, which was used as a control, has been published elsewhere (a gift from Prof. Tomomichi Fujita).

**RT-PCR** To test the expression of the *cals3m* transgene, total RNA was isolated using the PureLink RNA Mini Kit. RNA was treated with the ezDNase enzyme prior to reverse transcription with SuperScript IV VILO Master Mix, following the manufacturer's guidelines. PCR was subsequently performed with GoTaq G2 DNA polymerase using the primer pairs YCL-18F and YCL-18R, YCL-19F and YCL-19R, and YCL-19F and YCL20R pairs, targeting either the beginning or the end of the *cals3m* transgene. Expression of *PpUBI* was used as a reference to control for gene expression levels between samples. Thermocycling conditions were: 94°C for 2 min; 94°C for 30 sec, 58°C for 30 sec, 72°C for 60 sec (25 cycles); 72°C for 5 min. Primer melting temperatures were between 58 and 60°C.

**Quantitative RT-PCR** Lateral branches were dissected out with thin tweezers and shoot apices were cut off with micro-scissors. Total RNA was isolated using the PureLink RNA Mini Kit. RNA was DNase treated with
ezDNase enzyme prior to reverse transcription with SuperScript IV VILO Master Mix following manufacturer's guidelines. qRT-PCR was performed using a ROCHE Fast Start Universal SYBR Green master mix and a ThermoFisher StepOnePlus Real-Time PCR System. Thermocycling conditions were: 95°C 10 min; 95°C 10 sec, 60°C 30 sec (40 cycles). Specificity of PCR amplification was checked by melt curve analysis at each run. *PpGHL17-1* (Pp3c17_13760.1) and *PpUBI* (Pp3c12_5760V3.1) genes were respectively amplified using the following primer pairs: PDR4.1F, PDR4.1R, YCL-12F-R1-UBI and YCL-12R-R1-UBI. Expression level of *PpUBI* was used as a reference to normalize gene expression level between conditions, and Ct values were converted into relative expression values using the comparative Ct method⁶⁷.

**β-estradiol treatment** One-week-old moss colonies grown on BCDAT were immersed in a solution of 1 µM β-estradiol (C₁₈H₂₄O₂) for 24 hours in normal culture conditions before RNA extraction or aniline blue staining. The β-estradiol working solution was prepared from a stock solution at a concentration of 10 mM by performing a 1:10,000 dilution in deionized water. The control solvent was a 1:10,000 dilution of DMSO in deionized water.

**GUS staining and imaging** Gametophores were isolated from colonies grown on BCDAT for 4 weeks and incubated at 37°C in a 100 mM phosphate buffer (pH 7) containing 10 mM Tris (pH 8), 1 mM EDTA (pH 8), 0.05% Triton X-100, 5 mM potassium ferricyanide, 5 mM potassium ferrocyanide, and 2 mM X-GlcA (5-Bromo-4-chloro-3-indolyl-β-D-glucuronic acid) for 60-105 minutes. Tissues were then de-stained in 70% ethanol and imaged with a Keyence VHX-900F digital microscope with a 5-50X or a 50-200X objective. Three to four independent transgenic lines were observed for each genetic construct, and the pictures shown are representative of the GUS signal detected in the majority of gametophores.

**Aniline blue staining and quantification** Aniline blue staining and quantification were performed as described in the detailed protocol by Muller *et al.*³⁴. In brief, a fresh solution of aniline blue fluorochrome was prepared by a 1:3 dilution in 0.1 M phosphate buffer (pH 8.5, 97 mM K₂HPO₄, 3 mM KH₂PO₄) from a stock solution at 0.1 mg/mL aniline blue dissolved in sterile deionized water and stored at 4°C in the dark. Plant tissues were harvested from the edge of a 10-14 day old colony grown under standard conditions and immersed in 100 µl of aniline blue working solution for 30 minutes in the dark. Plant tissues were then rinsed with deionized water and aniline blue fluorescence was imaged using a Leica SP8 microscope confocal microscope (40X water immersion lens, excitation from the Diode 405 nm laser line, and emission collected at 440-520 nm). Z-stacks were processed with the image analysis Fiji to quantify the intensity of the fluorescence signal as described in Muller *et al.*³⁴.

## HPTS dye diffusion assay in whole gametophores

**Plant growth** Moss colonies were initiated from protonemal cultures and cultivated on BCD medium (250 mg/L MgSO₄.7H₂O, 250 mg/L KH₂PO₄ (pH 6.5), 1010 mg/L KNO₃, 12.5 mg/L FeSO₄.7H₂O, and 0.001% Trace Element Solution (0.614 mg/L H₃BO₃, 0.055 mg/L AlK(SO₄)₂·12H₂O, 0.055 mg/L CuSO₄.5H₂O, 0.028 mg/L KBr, 0.028 mg/L LiCl, 0.389 mg/L MnCl₂·4H₂O, 0.055 mg/L CoCl₂·6H₂O, 0.055 mg/L ZnSO₄.7H₂O, 0.028 mg/L KI and 0.028 mg/L SnCl₂·2H₂O)), supplemented with 7 g/L agar, 5 g/L sucrose and CaCl₂ added to a 1 mM concentration after autoclaving. They were grown for 4 weeks at 23°C under a 16h light/8h dark cycle, at 50-150 µmol.m⁻².s⁻¹ in MLR-352 growth cabinets, and then transferred to the dark at 23°C for an additional 6-8 weeks.

**HPTS bead preparation** Sephadex G50 beads were suspended in Milli-Q water to 1.25% (w/v) (0.625 g in 50 mL). After 2 min, beads were centrifuged for 2 min at 400 × g and the supernatant was removed. A volume of 500 µL of the beads was transferred to a 2 mL Eppendorf tube and mixed with 500 µL of 10 mM HPTS (8-hydroxypyrene-1,3,6-trisulfonic acid trisodium salt). The HPTS-bead mixture was centrifuged for 1 min at 10,000
× g, and the supernatant was removed. A second round of centrifugation (2 min at 10,000 × g) was performed, and the supernatant was again discarded.

**HPTS bead application and fluorescence imaging** Single gametophores were sampled from dark-grown colonies and placed vertically on 10 x 10 cm plates containing agar-supplemented BCDAT medium. Gametophore tips were cut off using micro-scissors, and an HPTS-impregnated Sephadex bead was applied to the stump of each gametophore using fine tweezers. Shoots were incubated for 2 hours, after which the beads were removed and gametophores were mounted on a microscope slide and covered with a coverslip. HPTS fluorescence was imaged using a Leica SP8 confocal microscope (HC PL Fluotar 10×/0.3 1 WD 11 mm lens; excitation at 448 nm; emission collected at 480–560 nm). The distance from the cut site to the front of dye diffusion was measured using Fiji software.

## Quantification of gametophores branching and geometry

**Plant growth** For quantification of gametophore branching patterns, phytomer length and stem diameter, moss colonies were cultivated on BCDAT for 5-8 weeks in sterile Magenta GA-7 tissue culture vessels under standard growth conditions, as described previously¹³,¹⁵.

**Branching patterns** For each genotype analyzed, sixty gametophores were isolated from the rest of the colony and dissected from the base to the tip to count the phytomer number and plot branch position.

**Phytomer length** Thirty-six wild-type gametophores were isolated and dissected from base to tip to count the number of phytomers, excluding the apical elongation zone based on visual clues. Gametophores were then imaged with a Keyence digital microscope (VHX-900F) with a 5-50X objective. The region where phytomers were counted was directly measured using the built-in software. The average phytomer length was estimated by calculating the slope of a linear regression (using `scipy.optimize.curvefit`) of the gametophore portion length as a function of the associated phytomer number.

**Stem diameter** Sixty-nine wild-type individual gametophores were isolated and dissected from base to tip to count the number of phytomers. The gametophores were then imaged with a Keyence digital microscope (VHX-900F) with a 50-200X objective. Using the built-in software, three measurements of stem diameter were taken from the apical, middle and basal region along the apical-basal axis and for each gametophore.

**P. patens** gametophore feeding with ¹³C₆-IAA stable isotope Six-week-old gametophores (approximately 50 mg per replicate) were incubated in 2 ml of liquid BCD medium supplemented with 1µM ¹³C₆-IAA stable isotope (Cambridge Isotope Laboratories, Tewksbury, MA, USA). Gametophores tissues were collected in three replicates at specific timepoints (0, 0.5, 2, 6, 12, 24 hours; T0 samples were rinsed in the medium for 10 seconds). After incubation, the medium was removed, and the tissues were placed on a nylon net filter and washed with dH₂O to remove residual medium while draining under negative pressure. The samples were then rapidly re-weighed, flash frozen in liquid nitrogen, and stored at -80°C.

**LC/MS analysis** Samples (ca. 20 mg, fresh weight) were homogenized with 1.5 mm zirconium beads using a FastPrep-24 instrument (MP Biomedicals, CA, USA) with 100 µl of 1M HCOOH in water. After centrifugation at 30,000 x g and 4°C for 20 minutes and re-extraction, the combined supernatant was applied to an SPE Oasis HLB 10mg 96-well plate (Waters, Milford, MA, USA). The supernatant was pushed through the SPE plate using a Pressure+96 manifold (Biotage, Uppsala, Sweden). The 96-well SPE plate was washed three times with 100
µl of dH₂O. Samples were eluted with 100 µl 50% acetonitrile/water (v/v). An aliquot of the SPE eluate was injected into the LC/MS system. Metabolites were separated on a Kinetex EVO C18 column (2.6 µm, 150 x 2.1 mm, Phenomenex, Torrance, CA, USA). The mobile phases consisted of A) 5 mM ammonium acetate and 2 µM Medronic acid in water and B) 95/5 acetonitrile/water (v/v). The following gradient program was used: 5% B at 0 min, 7% B at 0.1 to 5 min, 10 to 35% at 5.1 to 12 min, 100% B at 13 to 14 min, and 5% B at 14.1 min. Analysis was performed on an LC/MS system consisting of a UHPLC 1290 Infinity II coupled to a 6495 Triple Quadrupole mass spectrometer. MS analysis was performed in MRM mode. The following MRM transitions were used: ¹³C₆-IAA (-) 180; 136; ¹³C₆-oxIAA (+) 198; 152. Data acquisition and processing were performed using Mass Hunter software B.08.

## Plasmodesmata imaging and quantification

**Tissue preparation** Seven week-old gametophores were fixed in a solution containing 2.5% glutaraldehyde, 2.5% formaldehyde, 0.18 M sucrose, 1% tannic acid in 0.05 M sodium cacodylate buffer (pH 6.5). The tissues were then post-fixed overnight in 1% osmium tetroxide, rinsed with water, dehydrated using a graded ethanol series over 4 hours, and embedded in epoxy resin. Embedded fragments were sectioned using an EM UC6 ultramicrotome with a JUMBO diamond knife. Semi-thin sections (250–500 nm) were collected at regular intervals on glass slides and stained with toluidine blue in order to locate regions of interest in the tissues. Ultra-thin sections (80 nm) were then collected on formvar-coated copper grids for transmission electron microscopy (TEM) and on silica wafers for scanning electron microscopy (SEM). For SEM, sections were first stained with 2% uranyl acetate for 15 min, washed three times with water, subsequently stained with lead citrate for 6 min, followed by three additional washes with water. Sections were then mounted onto stubs with double-coated adhesive carbon tabs and sputter-coated (Safematic compact coating unit) with a 2 nm carbon layer, transparent to electrons, to prevent electron accumulation in the tissues.

**TEM and SEM imaging** Gametophore sections were imaged at high resolution using transmission electron microscopy (TEM) (sections A-E, Figures S2 and S3) or scanning electron microscopy (SEM) (sections F-G, Figures S2 and S3). SEM streamlined both image acquisition and post-processing (mosaic reconstruction) compared to TEM. Sections for TEM were imaged with a JEOL JEM-1400 microscope, operating at 120kV, at 2nm pixel resolution, equipped with a RIO9 or a Orius camera (Ametek). Mosaic images were created using GMS3 and IMOD software. Sections for SEM were imaged with a Zeiss Crossbeam 550 microscope operating at 5kV and 500pA. Image acquisition was performed with a BSD detector, at 2-10 nm pixel resolution, driven by the Atlas5 software, with several tiles covering the entire region of interest to generate a mosaic image. Mosaic images from SEM were reconstructed using Atlas5.

**Plasmodesmata maps** Plasmodesmata were manually counted on TEM and SEM images using the ROI manager in Fiji. The counts were converted to densities by dividing the plasmodesmata number by the cell wall length and the section thickness (80 nm). To generate the cell templates for mapping plasmodesmata densities, cell walls were manually drawn from each microscopy image using inkscape, and the drawings were converted into binary images which were segmented using the Python library timagetk. From the segmented image, an incidence graph representing the cell template was generated using the Python library cellcomplex.

**Estimation of intercellular diffusivity in protonemal filaments** To quantify intercellular diffusion via plasmodesmata, photoconversion experiments were performed in linear protonemal filaments of Dendra2-expressing lines. Protonemal tissue was grown in 35 mm glass bottom dishes containing BCDAT medium supplemented with
0.5% (w/v) glucose. After a 2-day incubation period under continuous light, the tissue was overlaid with a cov-erglass and incubated for an additional week under unilateral red light⁶⁸. Then, 24 h prior to photoconversion, 5 µM β-estradiol or a control solution containing an equal amount of solvent (0.05% v/v DMSO) was added to the tissue. Photoconversion of Dendra2 in the 10th cell from the apex was performed on a Nikon Ti-eclipse microscope body equipped with a scanhead consisting of two galvano-scanner mirrors. Local photoconversion of Dendra2 was achieved by scanning the 405 nm laser light (generated by a 120 mW Omicron LuxX diode laser operating at 10% of its maximum output) for 150 passes over a pre-defined region of the imaging plane using a 20x 0.75NA Plan Apo VC objective. After a 15-min post-photoconversion waiting period³⁵, t=0 images were taken using a 10x 0.3NA Plan Fluor lens and a Yokogawa CSU-X1 spinning disk head. Samples were then returned to incubation conditions and imaged again 16 h post-photoconversion. Green fluorescence of Dendra2 was imaged using 491 nm excitation light, with emission light collected through a 527/60 bandpass filter. Red Dendra2 fluorescence was imaged using 561 nm excitation light, with emission light collected through a 595/50 bandpass filter. Image digitization was performed with a Photometrics Prime 95B sCMOS camera equipped with a 1.2x post-magnification lens. All components were operated by MetaMorph software.

Quantification and statistical analysis

All statistical tests correspond to Mann-Whitney U two-sided tests unless specified otherwise.

Modeling and implementation

**Tissue representation** Both 2D (real and idealized sections) and 3D (idealized gametophore sections) cellular templates were generated using the Python library `cellcomplex` (https://mosaic.gitlabpages.inria.fr/cellcomplex/)⁶⁹. The resulting cell complexes consist of a set of topological elements of dimension 0 to 3, where higher-dimension elements are defined by the set of their lower-dimensional boundary elements. For 2D tissues, the faces (elements of dimension 2) represent the cells, each cell wall is defined by one to several edges (elements of dimension 1) depending on its curvature, and each edge is defined by exactly two vertices (elements of dimension 0), whose coordinates determine the spatial organization of the entire cell complex. For 3D tissues, the elements of dimension 3 represent the cells, each cell wall is defined by one to several faces depending on its curvature, and each cell edge is defined by one to several graph edges. Quantitative and qualitative properties may then be assigned to any element of any dimension, such as the cell wall type (longitudinal or transverse), plasmodesmata densities, source cell identity, phyllid axil cell identity and auxin concentration.

**Auxin diffusion equation implementation and solving** To numerically solve the auxin diffusion equation, we used a finite difference scheme and formalized the system in matrix form, incorporating no-flux (Neumann) boundary conditions at the tissue boundaries and constant concentration (Dirichlet) conditions in auxin source cells. We defined a square matrix $\mathcal{M}$ with dimensions corresponding to the number $N$ of cells in the tissue, such that

$$X' = \mathcal{M}X \quad (3)$$

where $X$ is the state variable, which is the vector of size $N$ of auxin concentration in each cell of the tissue. $\mathcal{M}$ is defined $\forall (i, j) \in \{1, \dots, N\}^2$ by $M_{ij} = 0$ if $i \in S$, otherwise

$$M_{ij} = \begin{cases} \mathcal{P}\rho(d_{ij})A_{ij}^{cw}/V_i & \text{if } j \in \mathcal{N}_i \\ -\mathcal{P}/V_i \sum_{j \in \mathcal{N}_i} \rho(d_{ij})A_{ij}^{cw} - \delta & \text{if } i = j \\ 0 & \text{otherwise} \end{cases} \quad (4)$$
with $\mathcal{P}$ the plasmodesmata permeability ($\mu m^3 h^{-1}$), $\rho(d_{ij})$ the plasmodesmata density ($\mu m^{-2}$) defined depending on the distance of the cell wall to the closest auxin source ($\mu m$), $A_{ij}^{cw}$ the area of the cell wall ($\mu m^2$), $V_i$ the cell volume ($\mu m^3$), $\delta(h^{-1})$ the auxin decay rate, and $S$ the set of auxin source cells. The system equilibrium $\mathbf{X}^*$ belongs to the kernel of $\mathcal{M}$ such that $\mathcal{M}\mathbf{X}^* = \mathbf{0}_N$. The dimension of the kernel of a matrix is equal to the number of rows in the matrix that consist entirely of zeros. Here, $\mathcal{M}_i = 0_N \forall i \in S$. Thus, the eigenvalue 0 has a multiplicity equal to the size of $S$, which is the number of auxin source cells. Moreover, for each auxin source cell $i \in S$, there exists a unique eigenvector $\mathbf{v}_i$ associated to the null eigenvalue that is non null in $i$, i.e. $\mathbf{v}_{ii} \neq 0$. Thus, the unique solution that satisfies the Dirichlet condition in auxin source cells is defined by a linear combination of these vectors, such that

$$ \mathbf{X}^* = \sum_{i \in S} \mathbf{v}_i \frac{\alpha}{\mathbf{v}_{ii}} \qquad (5) $$

Solving of eigenvectors and eigenvalues was done using the Python implemented function `numpy.linalg.eig`.

## Building an idealized gametophore

To generate a 3D idealized gametophore, we first needed to capture the key properties of cell geometry from the 2D images that we acquired. For this purpose, we created idealized two-dimensional longitudinal and transverse sections.

**Transverse section** To create an idealized transverse section, we used the function `circle_voronoi_topomesh` from the `cellcomplex` library. This function builds a circular voronoi diagram, a representation commonly used to model cell tissues due to its geometrical and topological properties. The function takes as inputs the average cell number along the radius of the cell complex and the average cell diameter. We tested various parameter values and compared the resulting cell complexes to our transverse histological sections in the mature gametophore region. Key metrics for comparison included the total area of the cell complex, the total number of cells, the average cell area, and cell wall lengths. Based on these comparisons, we selected a radius of 7 cells and an average cell diameter of 20 $\mu m$ (Figures S4A-S4C).

**Longitudinal section** To generate a virtual 3D cylindrical gametophore, we extruded the idealized transverse section over a number of cell layers. This idealization does not capture the conical shape of the apex and the staggered arrangement of cells. To anticipate the effect of these simplifications, we created the longitudinal section using a rectangular grid. For this idealization, we estimated the number of cells in width and height, as well as cell width and length. To maintain consistency with the transverse section, we fixed the cell width at 20 $\mu m$. Cell lengths were measured from the cell complex of the biological longitudinal section by projecting all vertices of each cell onto the inertia axis of the cell complex and calculating the distance between the two most distant projected points. The cell geometries in the longitudinal section displayed an elongation pattern, with cells in the apical region measuring 15 $\mu m$ long and progressively increasing to a median value of 80 $\mu m$. We assumed that this pattern could be described by a step function going from 15$\mu m$ to 80$\mu m$ long cells. Our analyses showed that setting the four most apical cell layers to a length of 15$\mu m$, followed by nine cell layers with a geometric increase up to 80$\mu m$ long cells, accurately reproduced the observed cell elongation pattern. Finally, we set the template width to 11 cells to approximate the biological section width, and the number of cell layers to 14 to approximate the biological section length (Figures S4D-S4F).

**Approximating the distribution of plasmodesmata densities as a continuous function** Our goal was to capture the observed plasmodesmata distribution and model it within our idealized cellular framework. To this end, we defined a continuous function accounting for two assumptions (based on our results, Figure 2) : (i) the
distance to the auxin source is an adequate proxy for the local plasmodesmata density, and (ii) plasmodesmata dilution on longitudinal walls is proportional to the wall area. This function was denoted as h, wherein

$$ \rho(d_{ij}) = \begin{cases} h(d_{ij}) & \text{for transverse faces} \\ h(d_{ij})A_{trans}^{cw}/A_{ij}^{cw} & \text{for longitudinal faces} \end{cases} \qquad (6) $$

where $d_{ij}$ is the distance of the cell wall at the interface between cells i and j to the closest auxin source cell, and $A_{trans}^{cw}$ is the average transverse cell wall area.

To capture the sigmoidal shape of the plasmodesmata distribution observed, we expressed h as a decreasing Hill function of coefficient 2

$$ h(x; a, K, c) = c + a \frac{K^2}{x^2 + K^2} \qquad (7) $$

with c denoting the basal plasmodesmata density, and a and K defining the amplitude and size of plasmodesmata enrichment around the source, respectively.

Then, we optimized the parameterization of h to reproduce the auxin distribution at equilibrium estimated on the biological tissues. This approach is more relevant than directly fitting the plasmodesmata distribution for two reasons. First, there is significant variability in plasmodesmata densities between cell walls. Since diffusion is a smoothing process, comparing its output (the auxin distribution) rather than its input (the plasmodesmata densities) allows us to capture the essential properties of the plasmodesmata distribution. Second, our cellular framework is idealized and lacks certain properties of the biological tissue, such as staggered cells and a conical apex. These topological properties contribute to the effective diffusivity in the tissue. By optimizing the plasmodesmata distribution on the idealized sections to reproduce the auxin distribution simulated on the biological tissue, we can predict an idealized plasmodesmata distribution that compensates for these simplifications.

To accomplish this, we simulated auxin diffusion from meristematic cells on the biological cell templates, assuming a minimum density of 0.2 CS/µm² for all cell walls (i.e., the minimal non-null density calculated among all sections) to ensure complete tissue connectivity and comparable outputs in model simulations. We subsequently estimated the fitting curve of the resulting auxin distribution, and employed it as a reference. The optimization was performed by screening the parameter space of h for a in [0, 100] with a step of 1, for K in [0, 50] with a step of 1, and c in [0, 4] with a step of 0.05. For each parameter set, we estimated the auxin distribution at equilibrium and its corresponding fitting curve, which we then compared to the reference. More precisely, the cost function was defined as the integral of the absolute difference between the two fitting curves over the tissue length. Fitting curves were estimated by optimizing the parameterization of the function $x \mapsto c + (1-c)K^2/(x^2 + K^2)$, which captures the auxin distribution shape and ensures to have the auxin concentration at the source ($x=0$) equal to 1 (Figure S5).

**Optimization on the idealized longitudinal section** To estimate the plasmodesmata distribution at the gametophore apex, we performed the optimization on the idealized longitudinal section, defining the entire first layer of the cell grid as an auxin source. This was done using the simulation on the biological tissue with the apical cell set as an auxin source. Optimal parameters were (32, 7, 3.1). However, based on a sensitivity analysis to different tissue idealizations (not shown), we set ($a_a, K_a, c_a$) to (30, 6, 3), which only marginally affected the results.

**Optimization on the idealized transverse section** To estimate the plasmodesmata distribution around the auxin source of a newly initiated branch, we performed the optimization on the idealized transverse section considering a single epidermal cell as an auxin source, compared with the simulation on the biological section at the level of a newly initiated branch with the branch apical cell as the auxin source. However, in this case, we could only compare plasmodesmata densities on longitudinal walls. Therefore, accordingly with our dilution
hypothesis, we optimized the parameters of $h(x; a^*, K, c^*)$ where $X^* = X \times A_{trans}^{cw}/A_{long}^{cw}$, $X = a, K$, where $A_{long}^{cw}$ corresponded to the modeled longitudinal cell wall length (i.e., the cell length). The optimal parameter set was $(a^*, K, c^*) = (47, 8, 0.4)$, which we adapted to $(a_b, K_b, c_b) = (188, 8, 1.6)$ for 80µm long cells.

**Modeled plasmodesmata densities in the basal region of the gametophore** We assumed that there was no plasmodesmata enrichment at the base of the gametophore ($a = K = 0$), and that plasmodesmata densities were similar in the basal region compared to the middle region of the gametophore ($c = 1.6$). Dilution of plasmodesmata densities was modeled similarly to other regions.

**Estimation of the Apical Inhibition Zone proper**

In our hypothesis framework, due to the apolar nature of auxin diffusion, the apical inhibition zone (AIZ) is determined by inhibition arising from auxin sources at the apex and at the most apical branch. Thus, the observed AIZ corresponds to a random variable determined by the sum of the inhibition reach from the apex (the AIZ proper, denoted AIZ*) and from the latest branch (IBIZ) (Figure S6A). To extract the distribution of the apical inhibition zone proper (AIZ*), we formalized the problem as follows:

Let $Y$ be the AIZ random variable, $X$ the IBIZ random variable and $Y^*$ the AIZ* random variable. We hypothesized that there is a uniform probability of measuring the AIZ at any time between two consecutive branch initiation events. Under this assumption, the relationship between the three discrete random variables is given by

$$Y = Y^* + \mathcal{U}(0, X) \qquad (8)$$

where $\mathcal{U}(0, X)$ denotes the discrete uniform probability distribution from 0 to $X$. Therefore, assuming $X$ and $Y$ to be independent, the probability distribution of $Y^*$ can be estimated by

$$\mathbb{P}(Y^* = y^*) = \sum_{y} \sum_{x} \sum_{i=0}^{x} \left( \frac{1}{x} \mathbb{P}(Y = y) \mathbb{P}(X = x) \right) \mathbb{1}_{y-i=y^*} \qquad (9)$$

This equation considers all possible combinations of the random variables ($X, Y$, and the uniform variable from 0 to $X$ for each value of $X$) and redistributes the frequencies according to the corresponding value of $Y^*$.

The experimental data provided empirical frequencies of apical and inter-branch inhibition zones that were used to estimate the probability distribution of their respective random variables. We found no evidence of a correlation between the apical and inter-branch inhibition zone sizes (Figure S6D), justifying the use of the above formula. In wild-type gametophores, the observed apical inhibition zones were mostly distributed between 15 and 19 phytomers, corresponding to a range of 4 phytomers. Since the IBIZ size distribution peaks at 4 phytomers, we expected the AIZ* to be centered around 15 phytomers, which was confirmed by the outcome of the above computation (Figures S6B-S6E). In this study, simulation outputs for the apical inhibition zone were systematically compared with the AIZ*, for both control and transgenic lines.

**Optimization of normal parameter distributions to reproduce biological data** To compare the model outputs with biological data, we assumed that $P/\delta$ and $\tau/\alpha$ varied stochastically according to normal distributions. We aimed to find the parameterization of the normal distribution of $P/\delta$ and $\tau/\alpha$ that would best reproduce the biological data. We explored mean and variance values for the two parameters. Each parameter set for the normal distribution of parameters corresponded to a specific distribution of inhibition zone sizes, estimated with 1000 draws. These distributions were then compared to the biological data using the Wasserstein distance, and the optimal parameter set was defined as the one minimizing this distance.
The Wasserstein distance is a metric defined between probability measures or distributions. Metaphorically, it considers each distribution as a pile of earth. The energy to transform one distribution into another depends on the amount of earth to be moved multiplied by the distance it has to be moved. The Wasserstein distance, also called the Earth Mover's Distance, estimates the minimum energy needed for this transformation.

The problem was formalized as follows: inputs are (1) the observed frequencies of an IZ size $\Phi_{obs} = (f_n)_{n \in \mathbb{N}}$ with $f_n$ the frequency of the IZ size $n$ and (2) the model prediction $F: \mathbf{x} \mapsto F(\mathbf{x}) \in \mathbb{N}$ which associates an IZ size with any parameter vector $\mathbf{x}$ in the parameter space $\Omega$. The aim of the approach is to find the multivariate normal distribution $N(\mathbf{x}; \mu, \Sigma)$ such that if $\Phi_{pred} = (\int_{\Omega_n} N(\mathbf{x}; \mu, \Sigma) d\mathbf{x})_{n \in \mathbb{N}}$ with $\Omega_n = \{\mathbf{x} \in \Omega | F(\mathbf{x}) = n\}$, then $d_W(\Phi_{pred}, \Phi_{obs})$ is minimized, with $d_W$ the Wasserstein distance. Here, $\mathbf{x} = (\mathcal{P}/\delta, \tau/\alpha)$, such that $\mu = \begin{pmatrix} \mu_{\mathcal{P}/\delta} \\ \mu_{\tau/\alpha} \end{pmatrix}$ and $\Sigma = \begin{pmatrix} \sigma_{\mathcal{P}/\delta} & 0 \\ 0 & \sigma_{\tau/\alpha} \end{pmatrix}$.

Optimizations were performed by screening the four-dimensional space defined by $(\mu_{\mathcal{P}/\delta}, \mu_{\tau/\alpha}, \sigma_{\mathcal{P}/\delta}, \sigma_{\tau/\alpha})$. Values explored for $\mu_{\mathcal{P}/\delta}$ spanned over [300, 4000] with an approximate discretization step of 85. Values explored for $\mu_{\tau/\alpha}$ spanned over [0.01, 0.99] with an approximate discretization step of 0.005. To ensure that the parameter values remained within realistic bounds, the standard deviations for each parameter were defined in relation to their respective means, such that only outliers of the corresponding normal distribution would fall outside the bounds of $\Omega = [300, 4000] \times [0.01, 0.99]$. Outliers were identified using John Tukey's test, defined as values outside the range $[Q_1 - k(Q_3 - Q_1), Q_3 + k(Q_3 - Q_1)]$, where $k = 1.5$ and $Q_1$ and $Q_3$ are the first and third quartiles of the distribution. With this constraint, the values explored for $\sigma_{\mathcal{P}/\delta}$ spanned from 20 to a maximum dependent on $\mu_{\mathcal{P}/\delta}$, with a discretization step of 20. The values explored for $\sigma_{\tau/\alpha}$ spanned from 0.005 to a maximum dependent on $\mu_{\tau/\alpha}$, with a discretization step of 0.005. Predicted IZ size distributions ($\Phi_{pred}$) were estimated from 1000 samples drawn from the corresponding normal distribution using the Python implemented function `scipy.stats.multivariate_normal`. Wasserstein distances were computed using the function `scipy.stats.wasserstein_distance`.

### Estimation of Dendra2 intercellular diffusivity in protonemal filaments
To quantify the cell-to-cell spread of photoconverted Dendra2 in protonemal filaments, the average intensity of red-fluorescing Dendra2 was measured at the moment of photoconversion (t=0h) and at t=16h in regions of interest (ROIs) covering the 8th to 12th cells from the filament tip, with the 10th cell being the 'donor' cell where the photoconversion was induced. The background signal intensity, measured immediately adjacent to each cell, was subtracted. Next, for each photoconversion and time point, fluorescence intensities were normalized to the total fluorescence intensity in all five measured cells. To compare the results, we used these signal intensities to predict the effective diffusivity for each acquisition. For this, we optimized the diffusion coefficient in a simple diffusion model along a linear file of cells to reproduce the observations at t=16h, using the observations at t=0 as the initial condition. The problem was formalized as follows:

For each sample, given $y_i(t)$ the signal intensity measured at cells $i = i_0 - 2, \dots, i_0 + 2$ at times $t = 0, 16$, normalized at each timepoint by the sum of the intensities in all cells, with $i_0$ the cell where Dendra2 photoconversion is induced, find $D$ such that if $\mathbf{x}(t)$ is a solution of

$$ \begin{cases} \frac{dx_i(t)}{dt} = D(x_{i-1}(t) - 2x_i(t) + x_{i+1}(t)) & \text{for } i = 1, \dots, N \\ x_i(0) = y_i(0) & \text{for } i = i_0 - 2, \dots, i_0 + 2 \text{ else } 0 \end{cases} \tag{10} $$

with $N = 25$, $i_0 = 10$ and no-flux boundary conditions (Neumann homogeneous), then $\sum_{i=i_0-2}^{i_0+2} (x_i(t) - y_i(t))^2$, $t=16h$, is minimized.

The underlying hypotheses are: 1) the degradation rate of the signal intensity is constant across time and space,
allowing normalizing to isolate the sole effect of diffusion, and 2) all cells are of similar size within and between samples. Optimization was performed using the Python implemented function `scipy.optimize.basinhopping`. Mann-Whitney U tests were conducted to compare the resulting sets of predicted diffusivities across the different genotypes and treatments (Figures S9D-S9F).

## Legends for supplemental files

**Data S1. Model parameter values and oligonucleotide sequences. Related to Figures 1, 4 and 6 and STAR Methods.** (A) Model parameter values estimation. (B) Parameter values defining the normal distributions of $P/\delta$ and $\tau/\alpha$ optimized to capture biological data. Abbreviations: WT, wild-type; IZ, inhibition zone; AIZ*, apical inhibition zone proper; IBIZ, inter-branch inhibition zone; BIZ, basal inhibition zone. (C) Oligonucleotides used in this study.

## References

[1] C. F. Kratochwil and R. Mallarino (2023). Mechanisms Underlying the Formation and Evolution of Vertebrate Color Patterns. *Annual Review of Genetics* 57, pp. 135–156.

[2] D. Barthélémy and Y. Caraglio (2007). Plant Architecture: A Dynamic, Multilevel and Comprehensive Approach to Plant Form, Structure and Ontogeny. *Annals of Botany* 99, pp. 375–407.

[3] Y. Coudert (2017). The Evolution of Branching in Land Plants: Between Conservation and Diversity. In Evolutionary Developmental Biology: A Reference Guide. L. Nuno de la Rosa and G. Müller, Springer International Publishing, Cham, pp. 1–17.

[4] O. Leyser (2003). Regulation of shoot branching by auxin. *Trends in Plant Science* 8, pp. 541–545.

[5] B. J. Ferguson and C. A. Beveridge (2009). Roles for Auxin, Cytokinin, and Strigolactone in Regulating Shoot Branching. *Plant Physiology* 149, pp. 1929–1944.

[6] P. Prusinkiewicz, S. Crawford, R. S. Smith, K. Ljung, T. Bennett, V. Ongaro, and O. Leyser (2009). Control of bud activation by an auxin transport switch. *Proceedings of the National Academy of Sciences* 106, pp. 17431–17436.

[7] J. Bertheloot, F. Barbier, F. Boudon, M. D. Perez-Garcia, T. Péron, S. Citerne, E. Dun, C. Beveridge, C. Godin, and S. Sakr (2020). Sugar availability suppresses the auxin-induced strigolactone pathway to promote bud outgrowth. *New Phytologist* 225, pp. 866–879.

[8] K. V. Thimann and F. Skoog (1933). Studies on the Growth Hormone of Plants. *Proceedings of the National Academy of Sciences* 19, pp. 714–716.

[9] Y. Wang, J. Wang, B. Shi, T. Yu, J. Qi, E. M. Meyerowitz, and Y. Jiao (2014). The Stem Cell Niche in Leaf Axils Is Established by Auxin and Cytokinin in Arabidopsis. *The Plant Cell* 26, pp. 2055–2067.

[10] M. A. Domagalska and O. Leyser (2011). Signal integration in the control of shoot branching. *Nature Reviews Molecular Cell Biology* 12, pp. 211–221.

[11] S. Yadav, H. Kumar, M. Mahajan, S. K. Sahu, S. K. Singh, and R. K. Yadav (2023). Local auxin biosynthesis promotes shoot patterning and stem cell differentiation in Arabidopsis shoot apex. *Development* 150, dev202014.

[12] A. N. Stepanova, J. Robertson-Hoyt, J. Yun, L. M. Benavente, D.-Y. Xie, K. Dolezal, A. Schlereth, G. Jürgens, and J. M. Alonso (2008). TAA1-mediated auxin biosynthesis is essential for hormone crosstalk and plant development. *Cell* 133, pp. 177–191.
[13] M. Thelander, K. Landberg, A. Muller, G. Cloarec, N. Cunniffe, S. Huguet, L. Soubigou-Taconnat, V. Brunaud, and Y. Coudert (2022). Apical dominance control by TAR-YUC-mediated auxin biosynthesis is a deep homology of land plants. *Current Biology* 32, 3838–3846.e5.

[14] T. A. Bennett, M. M. Liu, T. Aoyama, N. M. Bierfreund, M. Braun, Y. Coudert, R. J. Dennis, D. O’Connor, X. Y. Wang, C. D. White, et al. (2014). Plasma Membrane-Targeted PIN Proteins Drive Shoot Development in a Moss. *Current Biology* 24, pp. 2776–2785.

[15] Y. Coudert, W. Palubicki, K. Ljung, O. Novak, O. Leyser, and C. J. Harrison (2015). Three ancient hormonal cues co-ordinate shoot branching in a moss. *eLife* 4, e06808.

[16] V. M. R. Spencer, L. Bentall, and C. J. Harrison (2023). Diverse branching forms regulated by a core auxin transport mechanism in plants. *Development* 150, dev201209.

[17] T. Fujita, H. Sakaguchi, Y. Hiwatashi, S. J. Wagstaff, M. Ito, H. Deguchi, T. Sato, and M. Hasebe (2008). Convergent evolution of shoots in land plants: lack of auxin polar transport in moss shoots. *Evolution & development* 10, pp. 176–186.

[18] A. Paterlini (2020). Uncharted routes: exploring the relevance of auxin movement via plasmodesmata. *Biology Open* 9, bio055541.

[19] S. W. Wu, R. Kumar, A. B. B. Iswanto, and J.-Y. Kim (2018). Callose balancing at plasmodesmata. *Journal of Experimental Botany* 6, pp. 326–15.

[20] D. Ušák, S. Haluška, and R. Pleskot (2023). Callose synthesis at the center point of plant development—An evolutionary insight. *Plant Physiology*.

[21] A. Vatén, J. Dettmer, S. Wu, Y.-D. Stierhof, S. Miyashima, S. R. Yadav, C. J. Roberts, A. Campilho, V. Bulone, R. Lichtenberger, et al. (2011). Callose biosynthesis regulates symplastic trafficking during root development. *Developmental cell* 21, pp. 1144–1155.

[22] Y. Liu, M. Xu, N. Liang, Y. Zheng, Q. Yu, and S. Wu (2017). Symplastic communication spatially directs local auxin biosynthesis to maintain root stem cell niche in Arabidopsis. *Proceedings of the National Academy of Sciences* 114, pp. 4005–4010.

[23] X. Han, T. K. Hyun, M. Zhang, R. Kumar, E.-j. Koh, B.-H. Kang, W. J. Lucas, and J.-Y. Kim (2014). Auxin-Callose-Mediated Plasmodesmal Gating Is Essential for Tropic Auxin Gradient Formation and Signaling. *Developmental Cell* 28, pp. 132–146.

[24] C. Gao, X. Liu, N. D. Storme, K. H. Jensen, Q. Xu, J. Yang, X. Liu, S. Chen, H. J. Martens, A. Schulz, et al. (2020). Directionality of Plasmodesmata-Mediated Transport in Arabidopsis Leaves Supports Auxin Channeling. *Current Biology*, pp. 1–20.

[25] N. M. Linh and E. Scarpella (2022). Leaf vein patterning is regulated by the aperture of plasmodesmata intercellular channels. *PLoS biology* 20, e3001781.

[26] N. L. Mellor, U. Voß, G. Janes, M. J. Bennett, D. M. Wells, and L. R. Band (2020). Auxin fluxes through plasmodesmata modify root-tip auxin distribution. *Development (Cambridge, England)* 147, dev181669–12.

[27] K. W. Rogers and A. F. Schier (2011). Morphogen Gradients: From Generation to Interpretation. *Annual review of cell and developmental biology* 27, pp. 377–407.

[28] S. R. Yu, M. Burkhardt, M. Nowak, J. Ries, Z. Petrášek, S. Scholpp, P. Schwille, and M. Brand (2009). Fgf8 Morphogen Gradient Forms by a Source-Sink Mechanism with Freely Diffusing Molecules. *Nature* 461, pp. 533–536.

[29] J. E. M. Dickmann, J. C. Rink, and F. Jülicher (2022). Long-range morphogen gradient formation by cell-to-cell signal propagation. *bioRxiv*, p. 2022.04.11.487794.
[30] V. A. Grieneisen, B. Scheres, P. Hogeweg, and A. F. M Marée (2012). Morphogengineering roots: comparing mechanisms of morphogen gradient formation. *BMC Systems Biology* 6, p. 37.

[31] A. Östin, M. Kowalyczk, R. P. Bhalerao, and G. Sandberg (1998). Metabolism of Indole-3-Acetic Acid in Arabidopsis. *Plant Physiology* 118, pp. 285-296.

[32] L. Z. Drábková, P. I. Dobrev, and V. Motyka (2015). Phytohormone Profiling across the Bryophytes. *PLoS ONE* 10, e0125411.

[33] N. Kamamoto, T. Tano, K. Fujimoto, and M. Shimamura (2021). Rotation angle of stem cell division plane controls spiral phyllotaxis in mosses. *Journal of plant research*.

[34] A. Muller, T. Fujita, and Y. Coudert (2022). Callose Detection and Quantification at Plasmodesmata in Bryophytes. In Plasmodesmata: Methods and Protocols. Y. Benitez-Alfonso and M. Heinlein, Springer US, New York, NY, pp. 177-187.

[35] T. Tomoi, Y. Coudert, and T. Fujita (2022). Tracking Intercellular Movement of Fluorescent Proteins in Bryophytes. In Plasmodesmata, Springer, pp. 321-332.

[36] S. Gombos, M. Miras, V. Howe, L. Xi, M. Pottier, N. S. K. Jasemi, M. Schladt, J. O. Ejike, U. Neumann, S. Hänsch, et al. (2023). A high-confidence Physcomitrium patens plasmodesmata proteome by iterative scoring and validation reveals diversification of cell wall proteins during evolution. *New Phytologist* 238, pp. 637-653.

[37] L. Wegner and K. Ehlers (2024). Plasmodesmata dynamics in bryophyte model organisms: secondary formation and developmental modifications of structure and function. *Planta* 260, p. 45.

[38] M. Arnoux-Courseaux and Y. Coudert (2023). Re-examining meristems through the lens of evo-devo. *Trends in Plant Science*.

[39] B. E. Juniper and P. W. Barlow (1969). The distribution of plasmodesmata in the root tip of maize. *Planta* 89, pp. 352-360.

[40] C. M. Duckett, K. J. Oparka, D. Prior, L. Dolan, and K. Roberts (1994). Dye-coupling in the root epidermis of Arabidopsis is progressively reduced during development. *Development* 120, pp. 3247-3255.

[41] T. Zhu, W. J. Lucas, and T. L. Rost (1998). Directional cell-to-cell communication in the Arabidopsis root apical meristem I. An ultrastructural and functional analysis. *Protoplasma* 203, pp. 35-47.

[42] T. Zhu, R. L. O'Quinn, W. J. Lucas, and T. L. Rost (1998). Directional cell-to-cell communication in the Arabidopsis root apical meristem II. Dynamics of plasmodesmatal formation. *Protoplasma* 204, pp. 84-93.

[43] C. van der Schoot and P. Rinne. "The symplasmic organization of the shoot apical meristem". English. In: *Plasmodesmata, structure, function, role in cell communication*. Ed. by A. van Bel and W. van Kesteren. Germany: Springer, 1999, pp. 225-242. ISBN: 9783540651697.

[44] L. G. Tilney, T. J. Cooke, P. S. Connelly, and M. S. Tilney (1990). The distribution of plasmodesmata and its relationship to morphogenesis in fern gametophytes. *Development* 110, pp. 1209-1221.

[45] I. Kim, K. Kobayashi, E. Cho, and P. C. Zambryski (2005). Subdomains for transport via plasmodesmata corresponding to the apical-basal axis are established during Arabidopsis embryogenesis. *Proceedings of the National Academy of Sciences of the United States of America* 102, pp. 11945-11950.

[46] N. Gerlitz, R. Gerum, N. Sauer, and R. Stadler (2018). Photoinducible DRONPA-s: a new tool for investigating cell-cell connectivity. *The Plant journal: for cell and molecular biology* 94, pp. 751-766.

[47] E. E. Deinum, B. M. Mulder, and Y. Benitez-Alfonso (2019). From plasmodesma geometry to effective symplasmic permeability through biophysical modelling. *eLife* 8, pp. 105-40.
[48] C. Faulkner, O. E. Akman, K. Bell, C. Jeffree, and K. Oparka (2008). Peeking into Pit Fields: A Multiple Twinning Model of Secondary Plasmodesmata Formation in Tobacco. *The Plant Cell* 20, pp. 1504–1518.

[49] E. M. Bayer and Y. Benitez-Alfonso (2024). Plasmodesmata: Channels Under Pressure. *Annual Review of Plant Biology* 75, pp. 291–317.

[50] W. J. Nicolas, M. S. Grison, S. Trépout, A. Gaston, M. Fouché, F. P. Cordelières, K. Oparka, J. Tilsner, L. Brocard, and E. M. Bayer (2017). Architecture and permeability of post-cytokinesis plasmodesmata lacking cytoplasmic sleeves. *Nature Plants*, pp. 1–11.

[51] A. B. Murray (2007). Reducing model complexity for explanation and prediction. *Geomorphology* 90, pp. 178–191.

[52] M. Kitagawa, T. Tomoi, T. Fukushima, Y. Sakata, M. Sato, K. Toyooka, T. Fujita, and H. Sakakibara (2019). Abscisic Acid Acts as a Regulator of Molecular Trafficking through Plasmodesmata in the Moss Physcomitrella patens. *Plant and Cell Physiology* 60, pp. 738–751.

[53] M. L. Brault, J. D. Petit, F. Immel, W. J. Nicolas, M. Glavier, L. Brocard, A. Gaston, M. Fouché, T. J. Hawkins, J.-M. Crowet, et al. (2019). Multiple C2 domains and transmembrane region proteins (MCTPs) tether membranes at plasmodesmata. *EMBO reports* 20, e47182.

[54] J. Pérez-Sancho, M. Smokvarska, G. Dubois, M. Glavier, S. Sritharan, T. S. Moraes, H. Moreau, V. Dietrich, M. P. Platre, A. Paterlini, et al. (2025). Plasmodesmata act as unconventional membrane contact sites regulating intercellular molecular exchange in plants. *Cell* 188, 958–977.e23.

[55] M. G. Johnston, A. Breakspear, S. Samwald, D. Zhang, D. Papp, C. Faulkner, and J. d. Keijzer (2023). Comparative phyloproteomics identifies conserved plasmodesmal proteins. *Journal of Experimental Botany* 74, pp. 1821–1835.

[56] A. F. Sankoh and T. M. Burch-Smith (2021). Plasmodesmata and hormones: pathways for plant development. *American Journal of Botany* 108, pp. 1580–1583.

[57] E. E. Tee and C. Faulkner (2024). Plasmodesmata and intercellular molecular traffic control. *New Phytologist* 243, pp. 32–47.

[58] M. Kitagawa and D. Jackson (2022). A Forward Genetic Approach to Identify Plasmodesmal Trafficking Regulators Based on Trichome Rescue. In *Plasmodesmata: Methods and Protocols*, Springer US, pp. 393–407.

[59] O. Wartlick, A. Kicheva, and M. Gonzalez-Gaitan (2009). Morphogen Gradient Formation. *Cold Spring Harbor perspectives in biology* 1, a001255–a001255.

[60] E. Véron, T. Vernoux, and Y. Coudert (2021). Phyllotaxis from a Single Apical Cell. *Trends in Plant Science* 26, pp. 124–131.

[61] C. J. Harrison (2017). Auxin transport in the evolution of branching forms. *New Phytologist* 215, pp. 545–551.

[62] Y. Coudert, N. E. Bell, C. Edelin, and C. J. Harrison (2017). Multiple innovations underpinned branching form diversification in mosses. *New Phytologist* 215, pp. 840–850.

[63] D. Lang, K. K. Ullrich, F. Murat, J. Fuchs, J. Jenkins, F. B. Haas, M. Piednoel, H. Gundlach, M. Van Bel, R. Meyberg, et al. (2018). The Physcomitrella patens chromosome-scale assembly reveals moss genome structure and evolution. *The Plant Journal* 93, pp. 515–533.

[64] M. Kubo, A. Imai, T. Nishiyama, M. Ishikawa, Y. Sato, T. Kurata, Y. Hiwatashi, R. Reski, and M. Hasebe (2013). System for Stable β-Estradiol-Inducible Gene Expression in the Moss Physcomitrella patens. *PLOS ONE* 8, e77356.
[65] P.-F. Perroud, D. J. Cove, R. S. Quatrano, and S. F. McDaniel (2011). An experimental method to facilitate the identification of hybrid sporophytes in the moss Physcomitrella patens using fluorescent tagged lines. *New Phytologist* 191, pp. 301–306.

[66] D. R. Mallett, M. Chang, X. Cheng, and M. Bezanilla (2019). Efficient and modular CRISPR-Cas9 vector system for Physcomitrella patens. *Plant Direct* 3, e00168.

[67] K. J. Livak and T. D. Schmittgen (2001). Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2CT Method. *Methods* 25, pp. 402–408.

[68] M. Kitagawa and T. Fujita (2013). Quantitative imaging of directional transport through plasmodesmata in moss protonemata via single-cell photoconversion of Dendra2. *Journal of Plant Research* 126, pp. 577–585.

[69] G. Cerutti, O. Ali, and C. Godin (2017). DRACO-STEM: An Automatic Tool to Generate High-Quality 3D Meshes of Shoot Apical Meristem Tissue at Cell Resolution. *Frontiers in Plant Science* 8.

[70] J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld, B. Schmid, et al. (2012). Fiji: an open-source platform for biological-image analysis. *Nature methods* 9, pp. 676–682.

[71] J. D. Hunter (2007). Matplotlib: A 2D Graphics Environment. *Computing in Science & Engineering* 9, pp. 90–95.

[72] C. R. Harris, K. J. Millman, S. J. Van Der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, et al. (2020). Array programming with NumPy. *nature* 585, pp. 357–362.

[73] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. *Nature methods* 17, pp. 261–272.

[74] M. Ishikawa, T. Murata, Y. Sato, T. Nishiyama, Y. Hiwatashi, A. Imai, M. Kimura, N. Sugimoto, A. Akita, Y. Oguri, et al. (2011). Physcomitrella cyclin-dependent kinase A links cell cycle reactivation to other cellular changes during reprogramming of leaf cells. *The Plant Cell* 23, pp. 2924–2938.
A

<figure><img src="image_3.png" /></figure>

B

<figure><img src="image_4.png" /></figure>

C

<figure><img src="image_5.png" /></figure>

D

<figure><img src="image_6.png" /></figure>
A

<figure><img src="image_7.png" /></figure>

B

<figure><img src="image_8.png" /></figure>

C

<figure><img src="image_9.png" /></figure>

D

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

PD density (μm⁻²)

E

<figure><img src="image_12.png" /></figure>

F

<figure><img src="image_13.png" /></figure>

G

<figure><img src="image_14.png" /></figure>

(Trans.)

(Long.)

H

<figure><img src="image_15.png" /></figure>

Distance to apical cell (μm)
A

<figure><img src="image_16.png" /></figure>

B

<figure><img src="image_17.png" /></figure>

C

<figure><img src="image_18.png" /></figure>

D

<figure><img src="image_19.png" /></figure>

E

<figure><img src="image_20.png" /></figure>

Distance to source (µm)
A

<figure><img src="image_21.png" /></figure>

B

<figure><img src="image_22.png" /></figure>

E

<figure><img src="image_23.png" /></figure>

C

<figure><img src="image_24.png" /></figure>

F

<figure><img src="image_25.png" /></figure>

D

<figure><img src="image_26.png" /></figure>

G

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

PD density (µm⁻²)

Auxin concentration

H

<figure><img src="image_30.png" /></figure>

I

<figure><img src="image_31.png" /></figure>

J

<figure><img src="image_32.png" /></figure>

Wasserstein distance

K

<figure><img src="image_33.png" /></figure>

L

<figure><img src="image_34.png" /></figure>

M

<figure><img src="image_35.png" /></figure>

Predicted IZ size (phytomers)

N

<figure><img src="image_36.png" /></figure>

O

<figure><img src="image_37.png" /></figure>

P

<figure><img src="image_38.png" /></figure>

Biological data
Model predictions
<figure><img src="image_39.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>E</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>D</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>F</figcaption></figure>

<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /><figcaption>Predicted IBZ size (phytomers)</figcaption></figure>
<caption>A</caption>

CALS4

20/23

<figure><img src="image_46.png" /><figcaption>B</figcaption></figure>

CALS4

<caption>C</caption>

CALS12

28/28

<figure><img src="image_47.png" /><figcaption>D</figcaption></figure>

CALS12

<figure><img src="image_48.png" /><figcaption>E</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>F</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>G</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>H</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>I</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>J</figcaption></figure>

<figure><img src="image_54.png" /></figure>

Predicted *IZ* size (phytomers)

<figure><img src="image_55.png" /></figure>
<figure><img src="image_56.png" /><figcaption>A</figcaption></figure>

IBIZ size (phytomers)

<figure><img src="image_57.png" /><figcaption>B</figcaption></figure>

IBIZ size (phytomers)

<figure><img src="image_58.png" /><figcaption>C</figcaption></figure>

phytomer number

<figure><img src="image_59.png" /><figcaption>E</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>H</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>G</figcaption></figure>

<figure><img src="image_62.png" /></figure>

IBIZ size (phytomer number)

<figure><img src="image_63.png" /><figcaption>J</figcaption></figure>

<figure><img src="image_64.png" /></figure>

IBIZ size (phytomer number)

<figure><img src="image_65.png" /><figcaption>I</figcaption></figure>

<figure><img src="image_66.png" /><figcaption>K</figcaption></figure>

<figure><img src="image_67.png" /><figcaption>L</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>M</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>O</figcaption></figure>

<figure><img src="image_70.png" /></figure>

Predicted IZ size
(phytomers)

<figure><img src="image_71.png" /><figcaption>P</figcaption></figure>

stem diameter ○ (cell number)

<figure><img src="image_72.png" /><figcaption>R</figcaption></figure>

<figure><img src="image_73.png" /></figure>

Model simulations

Biological data

IBIZ rank