Differential biosynthesis and cellular permeability explain longitudinal gibberellin gradients in growing roots

Significance Growth hormones are mobile chemicals that exert considerable influence over how multicellular organisms like animals and plants take on their shape and form. Of particular interest is the distribution of such hormones across cells and tissues. In plants, one of these hormones, gibberellin (GA), is known to regulate cell multiplication and cell expansion to increase the rate at which roots grow. In this work, biosensor measurements were combined with theoretical models to elucidate the biochemical mechanisms that direct GA distribution and how these patterns relate to root growth. Our detailed understanding of how GA distributions are controlled in roots should prove a valuable model for understanding the makings of the many other hormone distributions that influence how plants grow.

1 Simulating the root-tip growth dynamics using a cell-based model The growth zone of the Arabidopsis root tip comprises three spatially distinct developmental zones: the meristem, the elongation zone and the mature zone. We focused on mature plant roots in which the sizes of the meristem and elongation zone were stable (noting that for Arabidopsis plants, this stabilization occurs approximately 6 days after germination in our growth conditions). Thus, in the frame of reference moving with the start of the elongation zone, growth is quasi-steady and the spatial distributions of cell lengths and growth rates are independent of time. Within the growth zones, cell growth is anisotropic, and thus, we consider cells to have a time-dependent cell length but constant radius. In prescribing the cells' growth dynamics, we followed the results of Van der Weele et al 2003, and take the cells' relative elongation rate (RER) to be constant within each developmental zone. Using the RER values, the cell lengths evolve according to the following standard formula these equations were solved using matlab ode solver ode45 during the time intervals between successive division events. Cells in the meristem elongate slowly and divide at approximately regular time intervals. Based on the data in Beemster and Baskin 1998, cell division is stochastic; the time interval between the time at which 1 each cell is created the time when it next divides is chosen randomly from a normal distribution with mean T c hours and standard error T cse . For simplicity, we supposed that the mean cell division rate is constant throughout the meristem. On division, the daughter cells have half the length of the parent cell and a GA concentration equal to the parent cell. We track the positions of the cell centres. Cells move away from the root tip due to elongation of the more rootward cells and the division event only occurs provided the cell centre is within the meristem region. When a cell leaves the meristem and enters the elongation zone, it ceases division and undergoes rapid elongation (according to formula (1.1)). As detailed in the subsequent section, the ODE representing the GA dynamics in the elongation zone requires us to record the length of each cell i at the time it leaves the meristem, L mi . Simulations are started with two cells, which grow and divide until a large number of cells are produced and the dynamics have reached a quasisteady state in the frame of reference moving with the root tip. In the results presented, we simulated the model until 150 cells are produced -increasing this number to 300 makes no noticeable difference to the model predictions.
2 Specifying the growth dynamics using data In order to accurately simulate the growth dynamics in our experimental conditions, we collected data to parameterise the root growth kinematics for Col0, ga20ox-1,-2,-3 and ga2ox q. Analysing growth kinematics in this way has a long history; similar analysis can be found in Beemster and Baskin, 1998 and references therein. We first fit a simple piecewise linear function to measurements of cell lengths in terms of distance from the QC. Since we are assuming that the RER is constant within the meristem and EZ, and that the division rate is constant in the meristem, it is appropriate to assume that the average cell length is constant in the meristem, increases linearly with distance through the EZ, before becoming constant at the end of the elongation zone (see equations (2.5) and (2.6) below). This piecewise linear function can be characterised by four parameters: the length of the meristem, x m , the length of the EZ, x EZ , the average cell length in the meristem, L m , and the mature cell length, L mat . For each case, these four parameters were fitted to the cell-length data using Matlabs inbuilt fminsearch algorithm using multiple starting points. The fitted piecewise linear distributions showed reasonable agreement with the cell length measurements in Col0, ga20ox-1,-2,-3 and ga2ox q, see Fig S2D- We note that 1/c is often called the total cell production rate, i.e. the rate at which new cells are created by the meristem. We can then calculate the average cell division rate, D, by dividing the total cell production rate by the number of meristem cells (given by x m /L m ): By the law of exponential growth, the number of cells at time t equals N 0 exp(Dt) (where N 0 is the number of cells at time 0). Since the cell number doubles during the average cell cycle duration (i.e. time interval between successive cell division events), T c , we therefore obtain (see Beemster and Baskin, 1998).
Having calculated c and T c , we can then calculate the RER in each zone as follows. Integrating equation (1.1), provides a formula for the cell length at time t 2 , l(t 2 ), in terms of the cell length at a previous time t 1 , l(t 1 ), and the relative elongation rate in the meristem, RER m using this formula and supposing that during the time interval between successive cell divisions, T c , a cell has doubled in length (i.e. l(t 2 ) = 2l(t 1 )) we can calculate RER m via In the elongation zone, the RER can be calculated via Using the above formulae, we obtained the model growth parameters for each case, given in Table S1.
3 Simulating the GA dynamics in the cell-based model of the growing Arabidopsis root tip Following the method described in the Supplementary material of Band et al 2012, we now derive a system of ordinary differential equations (ODEs) that describe GA synthesis, degradation and dilution accounting for the subcellular structure.
Within each cell, GA is present in both protonated and anionic forms within the vacuole, nucleus and cytoplasm. We assume that within each of these compartments diffusion is sufficiently rapid that spatial variations within a compartment can be neglected, and thus denote the GA concentration in the cytoplasm by [GA] i , the nucleus by [GA] nuci and the vacuole by [GA] vaci . Transport between the nucleus and the cytoplasm is passive, whereas transport between the cytoplasm and the vacuole is controlled by the tonoplast; as is typical in hormone transport models, we assume that the tonoplast is impermeable to anionic GA, but that protonated GA can passively diffuse across it. Thus, the rate at which GA moves between the cytoplasm and vacuole depends on the difference between the cytoplasmic protonated GA concentration, given by [GA] i (1 + 10 pHcyt−pK ) −1 , and the vacuolar protonated GA concentration, given by [GA] vaci (1 + 10 pHvac−pK ) −1 , where pH cyt and pH vac denote the pH of the cytoplasm and vacuole respectively, and pK denotes GA's dissociation constant. GA synthesis and degradation are taken to occur in the cytoplasm since enzymes involved in degradation and the later steps of biosynthesis reside here (Olszewski et al, 2002). Thus, the GA concentrations within each subcellular compartment for cell i can be described by where V nuci (t), V cyti (t) and V vaci (t) are, respectively, the volume of the cell nucleus, cytoplasm and vacuole of cell i at time t, P nuc is the rate at which GA moves between the nucleus and cytoplasm, A ton is the tonoplast area, P ton is the permeability of the tonoplast to protonated GA, σ is the GA synthesis rate and β is the GA degradation rate. To represent the lower synthesis in the meristem and higher synthesis in the elongation zone, we suppose that the GA synthesis rate, σ depends on the distance from the QC, x, as detailed below. We suppose, for simplicity, that transport across the tonoplast is rapid (i.e. taking P ton A ton to be large), so that the cytoplasmic and vacuole concentrations of protonated gibberellin can be treated as equal (see equation (3.1c)), and therefore, letting κ denote the ratio between the gibberellin concentrations in the vacuole and cytoplasm, we have In addition, we assume that transport between the cytoplasm and nucleus is rapid. Thus, we take P nuc to be large, which suggests that Under the limits implied by the above assumptions, the model (3.1) then simplifies to equation for the cytoplasmic and nuclear gibberellin concentration (this being the sum of equations 3.1 a-c).

Meristem equations
Simulating ( where R denotes the cell radius (which it is appropriate to take to be constant). Therefore, in the meristem, where we have used equation (1.1). Here the first term of the right-hand-side represents synthesis and degradation, whereas the second term represents dilution as the cytoplasm and nucleus expand.

Elongation-zone equations
Within the elongation zone, the increase in cell volume is due primarily to an increase in vacuole volume. Thus, as cells elongate, gibberellin within the vacuole dilutes, leading to an influx of protonated gibberellin from the cytoplasm which leads to a reduction in gibberellin concentration within both the cytoplasm and nucleus. We take the simplest assumption and suppose that the cytoplasm and nuclear volumes are constant, and therefore where L mi denotes the length of the cell i when it passes the boundary between the meristem and elongation zone. We note that these elongation-zone assumptions are the same as those in our previous model of the GA dynamics in the elongation zone (Band et al 2012). Thus, in the elongation zone, equation where α denotes the ratio between the GA synthesis rates in the elongation zone to meristem and we have simplified the last term using equation (1.1). Here, the first term on the right-hand-side represents synthesis and degradation, whereas the second term represents dilution due to expansion of the vacuole.
Equations (3.6) and (3.8) are coupled to (1.1) and are solved using matlabs ode solver ode45 in the time intervals between the division events.

Specifying the spatially varying synthesis rate
To investigate our hypothesis that the higher GA levels exhibited in the elongation zone reflects an increase in GA synthesis rate, we specify a synthesis rate that increases with distance from the QC, x, of the form where constant σ QC represents the synthesis rate at the QC, α and ξ are constants, and n is a natural number. With this Hill function form, the synthesis rate varies smoothly between a low value and a high value, with a more step-like form for higher values of n, occurring at ξ microns from the QC.

Predicting the nlsGPS1 distribution
Once the model has predicted the GA distribution, we use the following relationship suggested in Rizza et al, 2017 to calculate the corresponding distribution of nlsGPS1 parameter values are given in Table S1.

Simulating the exogeneous GA experiments
To simulate the exogeneous GA experiments, we extended our governing equations to include an extra term representing the rate at which the exogeneous GA enters the cells. We suppose that the rate at which GA enters the root is proportional to the external GA concentration [GA] dose , with the constant of proportionality referred to as the permeability of the root to GA, P root . Thus, the governing equations (3.6) and (3.8) become for elongation zone cells.
We considered several cases for the form of the permeability, P root . Firstly we considered a constant permeability, see Fig 3D and Fig S8A. Secondly, we considered the permeability to depend on the distance from the QC, x; considering the form where P QC , ζ and λ are constants (see Fig 3E,F, Fig S8B). This form represents a permeability that varies smoothly between a low value close to the QC and a higher value for cells far from the QC.
6 Parameter considerations 6.1 Specifying the synthesis and degradation parameters via parameter surveys and comparison with sensor data The model parameters were specified from our own data or estimates in the literature where possible (as summarised in Table S1). Parameters related to the growth dynamics were estimated from our own growth data ( Figure S2, Table S1), parameters in the relationship between the GA concentration and sensor data were characterised in Rizza et al 2017, and parameters related to the pH and dissociation constant were available in the literature (Kramer, 2004;Kramer 2006 ; Table S1). However, the appropriate spatial distribution of the GA synthesis rate and the value of the GA degradation rate are unknown as these will depend on the in vivo activity levels of the metabolism enzymes. We investigated parameter space by performing a simple parameter survey, considering uniform distributions of the parameter values that characterise the GA synthesis rate distribution: (6.1) α ∈ [0.0002, 0.0004, · · · , 0.002], ξ ∈ [25, 50, · · · , 250] σ QC = [0.00001, 0.00002, · · · 0.0001], and calculating the difference between the predicted and observed sensor data via where Data(x i ) denotes the sensor data which comprises N measurements at distance x i from the QC (for i = 1, 2, N). The predicted sensor data at x i is denoted Model(x i ) and is calculated from the model prediction (using linear interpolation over space). The parameter survey was performed for two values of the Hill function co-efficient: n = 2 which enables us to consider a smoothly increasing synthesis rate and n = 10 which provides an example of a step-like increase in synthesis rate. Furthermore, we consider two values of the GA degradation rate (since the ga2oxq sensor data motivates the choice of the small degradation rate, as discussed in the main text). In each case, we selected the parameter set [α, ξ, σ QC ] which minimises f (6.2).
In the case of n = 2, β = 0.05, the σ QC parameter obtained was on the lower boundary of the investigated parameter space, therefore the parameter survey was rerun to investigate a space with smaller σ QC ∈ [0.000002, 0.000004, · · · , 0.00002]. The resulting parameter sets are give in Table S2.

Specifying the permeability parameters via parameter surveys and comparison with sensor data
In the simulations of the exogeneous GA, similar parameter surveys were performed. For the case of constant permeability, we estimated the single parameter, P root = P QC (setting λ = ζ = 0), whereas for the spatially varying permeability, we estimated the values of the three parameters in the permeability distribution (5.3), [ζ, λ, P QC ]. In each case, we considered the following uniform distributions for these parameters: ζ ∈ [1, 2, · · · , 10], λ ∈ [50, 100, · · · , 500], P QC = [0.05, 0.1, · · · 0.5]. (6.3) We selected the parameter set (given in Table S1) which minimises the difference between the model predictions and data 20 minutes after the application of exogeneous GA (by calculating f as specified in (6.2)). We also investigated a spatially varying permeability with a higher Hill co-efficient equal to 10 to represent a more step-like increase in permeability; however, estimating the parameters [ζ, λ, P QC ] and calculating the minimum f , (6.2), in each case, we concluded that there is better agreement between model predictions and data with the lower Hill co-efficient of 2 (as stated in (5.3)).

Assessing the influence of stochasticity on the model predictions
The model is stochastic, due to stochasticity in cell division times. The synthesis rate of a cell is specified based on the distribution given in using the position of the cell centre. Hence if this cell divides, the daughter cells will have different synthesis rates, based on the positions of their respective cell centres. Thus, the stochasticity in division events affects GA synthesis and hence the predicted GA distribution. This process only has a minor effect, however, and we find very little difference in the predicted GA distribution for different simulation replicates (see Figure S5E). The results presented throughout are an average from 20 replicates (where the predicted GA and nlsGPS1 distributions are interpolated to a uniform grid prior to averaging).

Assessing the influence of the growth parameters
Although the growth parameters were estimated from our own measurements or the literature, we investigate here how the choice of these growth parameters affect the model predictions.
Increasing the RER specified for the meristem or elongation zone increases the rate of GA dilution in the respective zone, as well as increasing the rate at which cells traverse the developmental zones. As one would expect, with a larger meristematic RER, we predict lower GA concentrations throughout the growth zone ( Figure S5A), whereas a larger elongation-zone RER results in lower predicted GA concentrations in the elongation zone ( Figure S5B). The proportion of meristematic cell volume that is cytoplasm, γ, has been estimated in dividing cells in Arabidopsis as approximately 0.7 (Willis et al, 2016) (although we note that this value was measured in the shoot apical meristem, and as far as we are aware this has not been measured in the Arabidopsis root meristem). Since key GA metabolism enzymes reside in the cytoplasm (Olszewski et al, 2002), our equations show that γ influences the effective magnitude of synthesis and degradation. The model predicts that a lower γ would result in a small decrease in GA concentration throughout the elongation zone ( Figure  S5C).
Another key parameter is the ratio between the GA concentration in the vacuole to the cytoplasm, κ. Here, we follow Band et al, 2012 in assuming that only protonated GA can cross the tonoplast into the vacuole, resulting in a vacuolar GA concentration that is approximately 1/15 of the cytoplasmic concentration. However, it is possible that GA transporters may be present and enable anionic GA to enter the vacuole, a possibility that is suggested by the recent observation of substantial fluorescence-tagged GA in the vacuole (Shani et al, 2013), suggesting larger values of κ may be possible. Increasing κ would increase the rate of dilution in the elongation zone where cell elongation is thought to be predominantly due to vacuolar expansion. Thus, as expected, with larger values of κ the model predicts lower GA concentrations in the elongation zone ( Figure S5D). We see that with κ greater than around 0.2, the model predicts a peak in the GA concentration in the centre of the elongation zone with the GA concentration then reducing as cells move towards the mature zone. We note that this peak is not observed in the nlsGPS1 data suggesting that should higher values of κ be shown to be appropriate in the future even higher elongation-zone synthesis would be required to counteract the effect of dilution.   Table S1 for further details). (G) Root elongation rate for wild-type and GA metabolic mutant backgrounds calculated from the root length data shown in panel C.  for wild-type, ga2oxq (reduced degradation rate, assumed to be 10% of the wild-type value) and ga20ox1, ga20ox2, ga20ox3 (reduced synthesis rate, assumed to be 10% of the wild-type value).    Predictions of nlsGPS1 emission ratios for wild-type versus ga2oxq (reduced degradation rate, assumed to be 10% of the wild-type value) with constant permeability. (B) Predictions of nlsGPS1 emission ratios for wild-type versus ga2ox q (reduced degradation rate, assumed to be 10% of the wild-type value) with gradually increasing permeability as shown in Fig 3E. (C) 3D images of nlsGPS1 emission ratios and YFP fluorescence of ga2ox q roots before and 20 minutes after 0.1 µM GA4. (D) nlsGPS1 emission ratios for individual nuclei of ga2ox q mutant in relation to distance in µm from the root tip, before and after treatment with GA4. Curves of best fit and 95% confidence intervals are computed in R using local polynomial regression (Loess) via ggplot, with smoothing parameter span=0.75. Scale bar: 30 µm.