Auxin fluxes through plasmodesmata modify root-tip auxin distribution

ABSTRACT Auxin is a key signal regulating plant growth and development. It is well established that auxin dynamics depend on the spatial distribution of efflux and influx carriers on the cell membranes. In this study, we employ a systems approach to characterise an alternative symplastic pathway for auxin mobilisation via plasmodesmata, which function as intercellular pores linking the cytoplasm of adjacent cells. To investigate the role of plasmodesmata in auxin patterning, we developed a multicellular model of the Arabidopsis root tip. We tested the model predictions using the DII-VENUS auxin response reporter, comparing the predicted and observed DII-VENUS distributions using genetic and chemical perturbations designed to affect both carrier-mediated and plasmodesmatal auxin fluxes. The model revealed that carrier-mediated transport alone cannot explain the experimentally determined auxin distribution in the root tip. In contrast, a composite model that incorporates both carrier-mediated and plasmodesmatal auxin fluxes re-capitulates the root-tip auxin distribution. We found that auxin fluxes through plasmodesmata enable auxin reflux and increase total root-tip auxin. We conclude that auxin fluxes through plasmodesmata modify the auxin distribution created by efflux and influx carriers.

• After Xuan et al. (2016), auxin biosynthesis is increased 10-fold in the two outer lateral root cap layers, rather than just in the quiescent centre and initials in the Band et al. (2014).
• Rather than a zero auxin boundary condition in the outer cell layers (epidermis, cortex and endodermis) at the shootward boundary as in Band et al. (2014), we use a zero flux condition between these cells and their rootward neighbours as in Mellor et al. (2016), by setting their steady state values to be equal.
For a more detailed description of the model and these changes, and a description of the addition of plasmodesmata, see below.

S1.1 Spatial structure
As in Band et al. (2014) (Supplemental methods, section 1) the tissues are based on a 2-dimensional cellular structure of a cross-secion of the Arabidopsis root (approximately 500µm shootwards from the root tip) obtained using confocal microscopy with cell walls stained with propidium iodide (see Materials and Methods). In each case the roots had been crossed with the DII-VENUS nuclear-located yellow-fluorescent-protein auxin-response reporter (Band et al. 2012;Brunoud et al. 2012), and the software package SurfaceProject (Band et al. (2014), supplemental methods, section 1.1) was used to extract a 2D plane from each image stack with all the nuclei brought into the plane of focus.
Once the 2D image was obtained the cell segmentation software CellSet (Pound et al. 2012) was used to extract the position of all cell walls and cells, and quantify DII-VENUS nuclear fluorescence for every cell. Cell types were input manually using CellSet, according the the template shown in Figure S1.
The tissue is set up using Python code based on the Openalea framework (Pradal et al. 2008), and defined as a set of points in space representing the vertices of the tissue, S1 Development: doi:10.1242/dev.181669: Supplementary information which are linked to wall objects, which in turn are linked to objects representing the cell compartments. Wall compartments may be shared between two cells in the tissue (which we define as 'inner' walls), or if associated with a cell at the edge of the tissue it is only associated with that one cell (and defined as an 'outer' wall). Cell membranes are represented by way of a directed graph between the cells. Dyson et al. (2014) recently obtained measurements for cell wall thickness in various cell types and wall orientations in the Arabidopsis root tip. We assume a constant cell wall thickness (λ = 0.14µm) approximately consistent with these measurements, so that the 2D area of a given wall is equal to the length of that wall multiplied by λ. The area of the small vertex compartments are approximated as λ 2 . Finally, the area of a given cell is calculated as the area of the polygon defined by the vertices bounding the walls adjacent to that cell. PIN1,2,3,4 and 7 efflux carriers and AUX1 and LAX influx carriers are positioned on cell membranes according to a set of rules based on cell type, position and membrane orientation, as described in Band et al. (2014), Supplemental Tables 2 and 3. One such carrier localisation used in this paper is shown in Figure S2.
In the model where ectopic PIN1 is expressed in the pin2 following the observations of Omelyanchuk et al. (2016) (Figures S7 and S9) we place PIN1 on the shootward facing membranes of the last 20 and 25 rootward epidermal and cortical cells respectively (on each side of the root).

S1.2.1 Carrier mediated flux
The model is based on the ordinary differential equations defined previously by Band et al. (2014), Supplemental methods, section 2. In this previous model auxin flux is defined from cell compartments to adjacent wall compartments (and vice versa), to represent movement of auxin across cell membranes, either via influx and efflux carriers or via passive diffusion. Constants determining the directionality of the carriers and to model the acid trapping of auxin within the cytoplasm are used as defined previously Band et al. (2014).
Five possible flux carrier-mediated components exist between each wall and cell compartment: passive (P IAAH ), PIN dependent (P P IN ), AUX1 dependent (P AU X1 ), LAX dependent (P LAX ), and a background carrier mediated flux (P back ). P IAAH and P back are ubiquitous, while P AU X1 and P LAX are both multiplied by one or zero depending on whether AUX1 and LAX are respectively present or absent for a given cell membrane.  Table 3, and Figure S2, with P P IN scaled by the total number of PINs present for a given membrane, so that e.g. a membrane having PIN1,3 and 7 has three times the flux capacity as a membrane with just PIN2. As in Band et al. (2014) we denote the flux from the k th apoplastic compartment between cells i and j to cell i as J ijk and set this S2 Development: doi:10.1242/dev.181669: Supplementary information flux as: where A 1 , A 2 , A 3 , B 1 , B 2 and B 3 are constants (described in Band and King (2012) and given in Table S2), P P IN , P AU X1 , P LAX , P IAAH and P back the permeabilities (see Table  S2), [Auxin] i is the auxin concentration in cell i, and [Auxin] aijk the auxin concentration in apoplastic compartment ijk.

S1.2.2 Apoplastic flux
As before, to simulate auxin diffusion within the apoplast, the flux between adjacent wall compartments is simulated by considering the flux from each wall to and from two small vertex compartments representing the two ends of a given wall compartment. Denoting the flux from vertex l to apoplast ijk as J ijkl we have the following: where D cw is the diffusion coefficient of auxin in the apoplast, S ijk the length of compartment ijk, and [Auxin] vl the auxin concentration in vertex l.

S1.2.3 Plasmodesmatal flux
In the new model presented in this paper, in addition to movement across the cell membranes and within the apoplast, we consider another flux, via intercellular plasmodesmata.
Since the plasmodesmata (where present) essentially link the cytoplasm of adjacent cells, we assume these fluxes are directly between the cells themselves, without being mediated by wall compartments. The plasmodesmatal flux from cell j to cell i (J p ij ) is therefore defined as : where P plas is the permeability per plasmodesmata (µm 3 s −1 )and d ij the plasmodesmatal density (µm −2 ) between cells i and j (d ij ≡ d ji ). The plasmodesmatal density varies by cell type and wall orientation and is specified using TEM data given by Zhu et al. (1998) and reproduced in Table S1. Rutschow et al. (2011) give an experimentally measured estimate of the plasmodesmata permeability in the Arabidopsis stele in the longitudinal direction of 8 µm sec −1 . Combining this value with the measured plasmodesmata density in the anticlinal walls of the stele of 9.92 µm −2 given by Zhu et al. (1998) gives an estimate for the permeability per plasmodesmata of P plas = 0.806 µm 3 s −1 .

S1.2.4 Production and degradation
We also include in each cell a constant, cell-type dependent auxin biosythesis rate (α HIGH i in the QC, initials and outer lateral root cap, α LOW i elsewhere, and referred to generically S3 Development: doi:10.1242/dev.181669: Supplementary information Table S1: Plasmodesmata densities (µm −2 ) used in the model, taken from Zhu et al. (1998).
Pericycle cell layer included in stele. See Figure S1 for cell type template. as α i ), and a uniform, linear degradation rate (β). The auxin biosynthesis rate is unifiorm expect for the QC and initials (as in Band et al. (2014)), and in the two outer lateral root cap layers (after Xuan et al. (2016)), where it is set to be 10-fold higher than in the remaining cell types.

S1.2.5 Combined equations
Combining the fluxes given by Equations (1)-(3) and the production and degradation terms we have the following set of ODEs: where R i is the 2-D area of cell i, C i denotes the set of cells adjacent to cell i, N ij denotes the number of apoplast compartments between cells i and j, V ijk denotes the pair of vertex compartments adjacent to apoplast compartment ijk and W l denotes the collection of apoplast compartments, ijk, adjacent to vertex l. Since the cells are defined by a set of ordered 2-D coordinates outlining an irregular polygon we can calculate their area using a simple triangulation algorithm.

S1.2.6 Boundary and initial conditions
For cells in the stele (not including the pericycle) at the shootward boundary we assume a constant fixed supply of auxin from the shoot so there is a fixed boundary condition of auxin ([Auxin] b = 1) in these cells. For the remaining cells at the shootward boundary (i.e. pericycle, endodermis, cortex and epidermis) we assume a zero gradient boundary condition, so that steady state auxin in these cells is equal to the value in the adjacent cell in the same cell layer; i.e. for a given outer boundary cell o with rootward neighbour n we have, at steady state:

S1.3 Numerical methods
Since in this paper we only consider steady-state values, and given the system is linear we can compute the steady state directly by setting the derivatives to zero and rearranging the system in the form: where (if n is the total number of cell, wall and vertex compartments) J is the n × n matrix representing all of the combined fluxes and degradation terms, [Auxin] is the n × 1 vector of auxin concentrations in every compartment, and r the n × 1 vector of production rates and boundary conditions. The resulting linear system is then solved using the sparse matrix solver spsolve from the Python package Numpy.
where i = 1, 2, · · · , N labels the cells, and the final term in (7) represents the change in auxin concentration due to fluxes across the cell membrane as defined above.
As described in detail in Band et al. (2012) , if we scale these equations and suppose that complex formation occurs rapidly (i.e. the rate constants k a , k d , l a , l d and l m are relatively large), we can reduce the network model to a single equation for the DII-VENUS concentration: where we define the parameters p 1 , p 2 , p 3 and p 4 as: and where the b subscript represents the steady-state value of a given variable at the stele shootward boundary (see section S1.2.6), and [TIR1] T is the total (conserved) concentration of TIR1/AFB receptors in each cell. We use the parameter values p 3 = 0.91, Given we are calculating the steady state auxin in every cell directly, we can set the derivative of equation (12) to zero, and rearrange to obtain the following steady-state relationship: S6 Development: doi:10.1242/dev.181669: Supplementary information where the * superscript represents the steady state values of Auxin and DII-VENUS in a given cell.

S2 Simulation cases
The default set of parameters are given in Table S2. The cases given in the paper where the parameters are altered or the model is otherwise perturbed are described in more detail below.

S2.1 Permeability per plasmodesmata
For the simulations without plasmodesmata shown in Figures 1 and S4a-d, the model is as described above, with P plas set to zero, while in the remaining figures P plas is set to our estimated value of 0.806 µm 3 s −1 unless otherwise stated. Notable exceptions to this are in Figures 4 and 5 where both the addition of 0.6µM H 2 O 2 to a wild type root and the addition of DEX to a DEX inducible gsl8 knockout mutant are simulated by a doubling of P plas to 1.612 µm 3 s −1 .

S2.2 Plasmodesmata density
Plasmodesmata density are set according to the values (in Table S1) obtained from Zhu et al. (1998), except in Figure 3 where the density is set uniformly to 0.83µm −2 (low plasmodesmatal density, measured value for periclinal walls between lateral root cap and epidermis), 5.42µm −2 (medium plasmodesmatal density, measured value for anticlinal epidermal walls) or 12.58µm −2 (high plasmodesmatal density, measured value for anticlinal endodermal walls).

S2.3 Transport mutants
The simulations of the pin2 and aux1 mutants are implemented by setting the level of each respective transporter to zero on all cell membranes, with all remaining model parameters unchanged. For the ectopic PIN1 in pin2 simulations (Fig. S6,S8), based on observations from Omelyanchuk et al. (2016) we add PIN1 to the shootward membranes of the 20 most rootward epidermal cells and the rootward membranes of the 25 most rootward cortical cells on either side of the root, while knocking out PIN2 entirely as before.

S2.4 Model evaluation against experimental data
To compare model predictions for DII-VENUS with measured fluorescense we normalise each value with the minimum value in each case, then plot the difference between model and data for each cell (Figures 1e,l, 2e,i, 4d, 5d, S4d,j S14d,j). To quantify this comparison (Figure 2j, S9) we take the mean absolute difference between cells in the model and cells in the data, i.e.: where C is the set of all cell compartments (with |C| denoting the number of cells), and model i and data i denoting model DII-VENUS and data DII-VENUS respectively, each normalised by their minimum values.

S2.5 Total auxin in tissue
The total auxin in the tissue (Auxin T ) as shown in Figure 2k for varying values of P plas in the different simulated genotypes is calculated as: where C is the set of all cell compartments, W is the set of all wall compartments, R i is the area of cell i, S j the length of wall j and λ the cell wall width.

S2.6 Single-file simulations
To assess the role of plasmodesmatal fluxes within the individual tissue layers, we simulated auxin transport through a single file of cells, with PIN efflux carriers located on the downstream cell membranes. The model incorporates passive diffusion of protonated auxin across cell membranes, PIN-mediated transport of anionic auxin across cell membranes and plasmodesmatal diffusion of auxin between adjacent cell cytoplasms. Labelling the cells by i = 1, · · · , N , we let c i (t) denote the auxin concentration in cell cytoplasm i and f i (t) denote the auxin concentration of the apoplast region neighbouring cell cytoplasms i and i + 1, at time t. The auxin fluxes across each cell membrane are then given by where J cf i denotes the flux from cell i to apoplast region i and J f ci denotes the flux from apoplast region i − 1 to cell i (see Band and King (2012) for the derivation of these flux terms). The auxin fluxes through plasmodesmata from cell i to cell i + 1 are given by where P plas denotes the permeability per plasmodesmata and d denotes the density of plasmodesmata (which is taken to be a constant in this single-layer model).
The auxin dynamics are then governed by the following system of coupled ordinary differential equations (ODEs): where l denotes the cell length and λ denotes the apoplast thickness. These ODEs, (17), are simulated under the assumption that the concentration in cells i = 1 and i = N are held fixed, c 1 (t) = 1, c N (t) = 0, and all other concentrations are initially zero, c i (0) = 0 for i = 2, 3, · · · , N , f i (0) = 0 for i = 1, 2, · · · , N . Figs. 3g, S12 show the simulation results for parameter values N = 100, l = 20 µm, λ = 0.5 µm, d = 1 µm −2 and the remaining parameters equal to those given in Table S2.           pin2 and (c) aux1 models with and without plasmodesmata (PD), for a range of values of the model permeability parameters P AU X1 , P LAX , P P IN and P back . P AU X1 omitted for the aux1 model as it is redundant in that case. In each case the remaining model parameters are as given in Table 2, Supplementary Modelling information. Development: doi:10.1242/dev.181669: Supplementary information (c) aux1

S25
(a) w.t. In each case the other PIN permeabilites are equal to P P IN as given in Table 2, Supplementary Modelling information. Development: doi:10.1242/dev.181669: Supplementary information (c) aux1  Figure S18: Simulated NPA treatment at four levels of efficacy in blocking PIN and background carrier efflux (25%, 50%, 75%, 100%, where %100 represents a complete block in efflux activity), in a model without plasmodesmata (A), and with plasmodesmata (B). The experimentally measured DII-VENUS following 2-hour treatment with 2 µM NPA is shown (on the same scale as the mock treated roots in Figure 1F by way of comparison), along with a normalised cell-by-cell comparison of the data versus the model with an estimated 75% reduction in carrier-mediated efflux. (C). Quantification of model fit versus percentage reduction in efflux carrier efficiency for models with and without plasmodesmata. Development: doi:10.1242/dev.181669: Supplementary information Auxin flux

S29
Auxin concentration DII-VENUS concentration P plas = 10 P plas = 4 P plas = 2 P plas = 1 Figure S19: Effect of increasing plasmodesmatal permeability (P plas ) on predicted auxin flux (top), auxin concentration (middle) and DII-VENUS (bottom). The colour scales on each set of images are set the same for easier comparison between values of P plas .  Figure S20: Horizontal cross sections of auxin concentration at various vertical positions in the root tissue for heterogenous plasmodesmata model (see Figure 2A for distribution), and models with uniformly low (0.83 µm −2 ), medium (5.42 µm −2 ) and high (12.58 µm −2 ) plasmodesmatal densities.  Figure S21: Root tip auxin distribution with uniform plasmodesmata density in pin2 (A-F) and aux1 (G-L). The predicted differences shown are between the auxin concentrations for the respective uniform plasmodesmata models and the variable plasmodesmata model (as shown in Figure 2F for pin2 and Figure S4G Figure S22: Effect of plasmodesmata on auxin propagation through a single file of cells. We suppose that auxin moves across cell membranes via both passive diffusion of protonated auxin and active transport mediated by PINs that, when present, are located polarly on the downstream membrane face of each cell. We suppose that auxin also passively diffuses between adjacent cell cytoplasms through plasmodesmata (in the case where plasmodesmata are present). See SI section 2.6 for the model equations. (P plas = 10 µm 3 s −1 ) DII-VENUS + 0.6 mM H 2 O 2 Figure S23: Replicates of background DII-VENUS distribution following 0.6 mM H 2 O 2 treatment. Scale bar 50 µm.  Figure S24: Effect of 0.6 mM H 2 O 2 treatment on the distributions of pin2 (a-f) and aux1 (gl) auxin and DII-VENUS. (a,g). Predicted steady-state auxin distribution (b,h). Predicted auxin fluxes (c,i). Predicted DII-VENUS distribution (d,j). Difference between predicted and observed DII-VENUS distribution (from predictions in panels c,i and data in panels e,k). (e,k). Quantification of DII-VENUS distribution using images in panels f and l. (quantified using CellSet image segmentation software). (f,l). Representative DII-VENUS confocal images. Scale bars 50 µm.  Development: doi:10.1242/dev.181669: Supplementary information