Rapid transporter regulation prevents substrate flow traffic jams in boron transport

In Arabidopsis thaliana, boron is shuttled from the soil across the epidermis, cortex and endodermis into the stele, to finally be taken up into the xylem, mediated by the action of secondary active boron exporters (BORs) and the enhancement of permeability through NIPs. Previously, we analysed boron patterning by using a two-dimensional cross-sectional model of the entire Arabidopsis root meristem, considering all spatial nuances of transporter localisation, polarity and intensity, whilst neglecting transporter dynamics and regulation (Shimotohno et al., 2015-04). Instead, transporter levels and activity were static, fixed to the observed steady state transporter expression under constant boron medium condition of 0.3 μM. While neglecting regulation of transporter expression or activity, a characteristic boron pattern emerged on the level of the root, presenting highest concentrations at the stem cell niche. This boron profile gradually decreased longitudinally up to the start of the differentiated tissues. The computationally simulated characteristic boron distribution, which we confirmed experimentally through LA-ICP-MS, strongly indicates that the mature root tissue is involved mainly in transporting boron from the soil into the xylem, while the distal meristematic tissues have a differential boron transport function, namely, to locally provide higher boron levels for fuelling local growth (i.e., cell wall material) at those regions. In that study, however, we did not address the dynamics of the transporter regulation nor how nutrient homeostasis is achieved. Here, we investigate the dynamics of the transporters and their mutual feedback with the generated boron distribution. We do so by spatially focussing on the differentiated tissue region involved in xylem-uptake.

At the elongation and differentiation zone of the root, transporters possess a striking polar expression; BORs are located at the inside-facing membranes of the cells, while NIPs are concentrated to the outside facets (see Figure 1A,B). Although protein levels can present strong variation between cell files, the transcription of the transporters takes place throughout the entire root tissue, even for NIP5;1 (Figure 1—figure supplement 1). We therefore use the simplifying assumption that all cell files are intrinsically the same in respect to the potential of expressing the transporters, and variations arise solely as a consequence of the nutrient distributions.

Figure 1 with 2 supplements see all

Download asset Open asset

Our analysis focuses on the transversal nutrient flow through the root that results from boron entering and leaving the different cell files transversally. Effectively, in our simple model, we only consider a single row of cells and, for simplicity, only four cell files over which boron is transported from the soil into the xylem (Figure 1C). To capture the dynamics of boron transport and transporter activities in root tissue, the model’s variables are the boron concentration in cells (C) and cell walls (W), and transporter activities of NIPs (N) and BORs (B) for each individual cell (n). The mutual dependency between these variables is described using ordinary differential equations (ODEs). For each cell file of the model, the transporters are produced at an equal rate, have the same transport potential, and are regulated in the same way. The modelled cells loosely map to the outermost epidermis (cell 1), the cortex (cell 2), the endodermis (cell 3), and finally to the pericycle/vasculature (cell 4). The last cell of the row is endowed with an upward convective flow, which captures the shootwards convective flow of the xylem. For reasons of simplicity, we consider these four cells to have equal dimensions, and make the same assumption for all intermediate cell walls (Figure 1C).

Given that transporter regulation can be observed in response to varying the boron conditions in the medium (Takano et al., 2005, Takano et al., 2006-06Takano et al., 2006-06, 2010; Miwa et al., 2013), (Figure 1—figure supplement 2), boron sensing underpinning those responses could either involve measurements of the intracellular concentration or involve measurements of the intercellular concentration (i.e., the levels within the apoplast). In the case of NIP5;1, we recently revealed that NIP5;1 mRNA degradation is triggered by boron-dependent ribosome stalling during the translation process (Tanaka et al., 2016). The boron-dependent degradation was also observed in an in vitro system without cell wall fraction, which suggests that the sensing mechanism depends on the cytoplasmic boron concentration. For the BORs the precise subcellular location of the boron sensing has not yet been determined, but it is reasonable to assume that the regulation is also responding to cytosolic boron levels, given that the purpose of the nutrient homeostasis is to keep the nutrient levels within the cell within bounds, and to this end directly sensing the cytoplasmic boron concentration is beneficial. In short, given the current knowledge, cytosolic boron sensing is the most parsimonious assumption.

Hence, our model assumes that all the boron transporters of a cell are regulated by the cytosolic boron concentration within that respective cell. Note that we loosely use the term ‘boron’ in our model to refer to both boric acid as well as borate. They can be found in a dynamic chemical equilibrium, ${\displaystyle B(OH{)}_3+H_{2}O{\displaystyle ⟺}B(OH{)}_4^{-}+H^{+}}$, with a ${\displaystyle pK}$ value of 9.24 (Woods, 1996). Given that this process does not involve any enzyme kinetics, this distinction is not explicitly considered in our model. Moreover, we solely focus on the import of soluble boron (both boric acid and borate), not considering any boron that is bound to the cell wall.

The temporal dynamics of the soluble boron concentration within the cytosol (${\displaystyle C_i}$) and in the cell walls (${\displaystyle W_i}$, being the cell wall adjacent to the inward facing facet of cell $i$) is determined by their inflow and outflow from and to neighbouring compartments:

Here, $N_i$ and $B_i$ represent boron permeabilities due to, respectively, the bi-directional transport by NIPs and the directional transport by BORs within cell $i$. Transporter-independent boron permeability through the plasma-membrane is also taken into account, described by parameter $p$. To capture the tissue context, boundary conditions are set at both extremities of the transversal cell series. The outer cell wall of the first cell (cell 1) is in contact with medium. Given the rapid diffusion rate of boric acid in water, we assume that the boron concentration in the outermost cell wall ($W_0$) is the same as in the medium, constant over time. The xylem transport that occurs from roots to upper tissue is represented by a removal term, with rate $a$, attributed to the innermost cell ($i=n$) of the cell row:

Assuming that transport activities are proportional to protein concentration, it follows that the dynamical regulation of NIP and BOR permeability ($N_i$ and $B_i$) are in direct proportion to the dynamics of their proteins (we will later show that our main insights do not depend on the assumption that the BOR transporter does not saturate under the range of concentrations here treated). It was therefore not necessary to introduce explicit variables for protein concentrations. Instead we directly use permeability, representing the protein concentrations and their regulation. The protein concentration, and thus the effective permeability, is determined by production and degradation rates.

As mentioned, accumulation of NIP5;1 mRNA is regulated through boron-dependent mRNA degradation. This suggest that the production of NIP5;1 protein is boron dependent. On the other hand, degradation of NIP5;1 protein is not found to be boron dependent (Takano et al., 2010; Tanaka et al., 2016). Thus, the levels of $N_i$ vary due to a production term that is dependent on cytosolic boron concentration and a constant degradation term:

where ${\alpha }_N$ is the production rate, ${\xi }_N$ is degradation rate, and $n_N$ is the Hill coefficient capturing the boron-dependent inhibitory regulation.

Accordingly, given that BOR protein levels are regulated through boron-dependent degradation (Takano et al., 2005), but there is no evidence of boron dependent production, $B_i$ varies over time due to constant production and boron-dependent degradation terms:

where ${\alpha }_B$ is the production rate, ${\xi }_B$ the basal degradation rate, $d_B$ the maximum boron-independent degradation rate, and $n_B$ the Hill coefficient capturing the boron-dependent inhibitory regulation.

We wish to use this model to understand the significance, if any, of the swift transporter regulation. To assess the effects of regulation swiftness on transport, we therefore introduce a time-scaling factor $ϵ$, which multiplies both the entire $\dot{N_i}$ and $\dot{B_i}$ terms (Equations 6 and 7). To align ourselves with the current experimental evidence that indicates swift regulation of NIP5;1 on the mRNA level and of BORs on the protein level, we modulate the production rate and degradation rate for NIP5;1 and BOR1, respectively. However, changing production rates of NIPs through a factor ${\displaystyle ϵ}$ without changing the degradation rate would result in changes of the overall protein concentration which might influence the system’s behaviour, not due to the swiftness of the regulation but due to the changes in equilibrium protein concentration. This also holds for the BOR regulation. Thus, to avoid such possible undesired effects, the time-scaling factor multiplies the entire right-hand side of the dynamics:

We define $ϵ=1$ as representing the swiftness of regulation in the wild type scenario, while reducing $ϵ$ allows us to simulate behaviours that would occur with less rapid regulatory dynamics, and increasing $ϵ$ allows us to consider even higher swiftness of regulation than experimentally observed. Note that the regulation of the transporter activity by cytosolic boron concentrations occurs through a sigmoidal relationship, in which $k_B$ and $k_N$ are the boron concentrations at the half maximum of the Hill function, defining the sensitivity of the response to the boron concentration. We set ${\displaystyle k_B\simeq 30k_N}$, consistent with our observations that the sensitivity of BOR expression is much less than that of NIP5:1 (Figure 1–Figure supplement 2), but will also show that our results do not require large differences between these two parameters.

We have assessed what is the biologically acceptable range for each parameter, based on the available literature (see Material and methods for details). Default parameter values were chosen to lie within these valid ranges, and are given in Table 1.

When regulation swiftness lies within the experimentally observed regime (${ϵ}_B={ϵ}_N=ϵ=1$), a constant flow of boron that is transported across the root ensues (Figure 2). Boron concentrations within the cells ($C_n$) remain constant over time. This steady state behaviour, however, is disrupted when the transporter regulation swiftness is decreased ten-fold ($ϵ=0.1$). Oscillations in the boron concentration arise in all cells, but most strikingly in the endodermis (cell 3) and cortex (cell 2) (Figure 2A), as well as in the cell walls (Figure 2B). Decreasing the transporters’ regulation swiftness even further ($ϵ=0.01,\mathrm{\hspace{0.17em}0.001}$) consistently enhances the amplitude of such oscillations, and enlarges their periods (Figure 2A,B). We conclude that even when environmental boron levels ($C_0$) and xylem activity ($a$) remain constant, simply deviating from the rapid experimentally observed transporter’s regulation swiftness is sufficient to push the boron transport system into an oscillatory dynamics.

Figure 2 with 1 supplement see all

Download asset Open asset

The reason the cells undergo such dramatic concentration variations at lower regulation swiftness can be intuitively linked to the accompanying changes in transporter levels which also ensue (Figure 2C,D). They are initially triggered by minute fluctuations in the concentrations, but then cause increasingly larger changes in cytosolic and apoplastic concentrations. It is not immediately clear, however, how it depends on either the NIP or BOR transporter regulation, why regulation swiftness exacerbates oscillation amplitudes and how their interdependencies either produce stable flows or unstable oscillations that propagate throughout the tissue. What we can already conclude though, is that within a polarised tissue rapid transport dynamics are necessary to ensure flux homeostasis and avoid instabilities in concentration levels. But to fully understand the process we will first observe how in the model the steady flow regime behaves in terms of transporter expression levels when we assume that the regulatory dynamics are as swift as experimentally observed.

Given that our question is focused on the possible advantages of swift transporter regulation within a polarised tissue context, our model purposely ignored details of tissue-specific regulation, using the simplifying assumption that all cells are endowed with the same potential of expressing transporters. This assumption may initially seem at odds with the actual distinct GFP levels of the GFP-tagged proteins (Figure 1B), which qualitatively indicates stronger BOR1 expression in the epidermis and endodermal files and pronounced NIP5;1 expression in the outermost cell file (lateral root cap or epidermis) only. It was therefore interesting to observe that our parsimonious model could already account for general features of the the observed transporters’ expression patterns (Figure 2C,D) solely as a consequence of the relative positioning of the cells within the cell row: in the model, NIP only becomes highly expressed in the outermost cell, the epidermis ($N_1$, representing NIP5;1 in cell 1, is by far the highest variable in the steady flow case). These relative levels are in qualitative agreement with the experimental findings regarding NIP5;1. Moreover, in contrast to this strong expression in the epidermis, in the model NIPs are present at much lower amounts in all the other inner cell files (cell 2, cell 3 and cell 4), again qualitatively similar to the experimental observations. These model results can be understood as follows. Although all in silico cells share the same intrinsic properties, as a consequence of NIPs acting as permeability facilitators rather than directional transporters, it is impossible for the cell file in contact with the medium (cell 1, the outermost cell file, that is, epidermal cell) to reach higher concentrations than the level in the medium. Only the next cell inwards, cell 2, is capable of reaching higher levels, due to the directed polar action of the BORs. The inevitable low boron levels in the epidermis give rise to the observed lack of NIP downregulation within the first cell file, and hence to its preferential expression.

The expression patterns of the active exporters, BORs, are quite different from NIP5;1, but again there are similarities between the model results (B values) and confocal images, with BORs distinctly expressed in the epidermis (cell 1) at the highest levels, in the cortex (cell 2) at high levels, whilst in the endodermis it is relatively weaker (cell 3) (Figure 1B). The model’s qualitative correspondence between the transporter expression patterns and those experimentally observed suggests that concentration-dependent regulation is at least sufficient for such expression patterns to arise. Corroborating this, GUS expression under the NIP5;1 promoter is indeed present in the inner root tissues, albeit the NIP5;1 protein levels in those cells are very low, supporting the notion that these files are likely responsive and capable of expressing NIP5;1, but inhibited to do so by the boron levels that arise due to the cells’ positioning within the tissue (Figure 1—figure supplement 1). However, these results, while suggesting sufficiency of our parsimonious assumptions in relation to observed transporter levels, do not offer evidence that tissue-specific regulation is not occurring. At most, the qualitative matches justify, at first instance, that the model is kept simple, without the need of introducing ad hoc tissue-specific dependencies, when exploring the effects transporter regulation dynamics. Furthermore, these results indicate that under our experimentally observed parameter values, the system presents steady nutrient throughput and constant cellular concentrations and transporter levels over time.

We next investigated more carefully the spatial-temporal characteristics of the oscillatory behaviour that arises when transporter swiftness is decreased (Figure 3A). Kymographs make it visually clear that oscillations in the individual cells are in fact spatially coupled, manifested as a boron wave that propagates backwards over the tissue, that is, the pulse moves in the direction contrary to the nutrient flow itself. The initiation of the instability, as seen by the first peak in concentration, occurs close to the xylem (i.e., the last cell, here, cell 4), followed by a strong rise in endodermal concentrations (in cell 3) around 90 min later. Due to boron’s inhibitory influence on the transporters, this intracellular rise leads to a decrease in both NIPs and BORs within that cell. This transporter downregulation causes a drop in boron throughput across that cell, leading to an accumulation of boron in the adjacent, outward facing, cell wall. Background permeability rates along the plasma membrane allow this apoplastic rise in concentrations to trigger an increase boron concentrations in the next outermost neighbouring cell, in this case, the cortex cell (cell 2). Cortex concentration levels thus rise, again triggering a shutdown of transporters, causing the same process to ensue in a spatially coupled manner in the outermost cell, the epidermis. In the epidermis, however, boron levels can never rise above the soil concentration levels, as discussed previously, such that the backward travelling wave loses its amplitude and terminates. The overall process is observed to be cyclical, with levels rising again close to the xylem, until a next wave is triggered in the innermost cells (Video 2). Moreover, we found that such traffic-jam-like behaviour is robust when smaller differences between $k_B$ and $k_N$ are considered (Figure 2—figure supplement 1A–C), as long as $k_B$ is sufficiently large as not to too strongly interrupt boron throughput through the tissue altogether (Figure 2—figure supplement 1D). Obviously, if $k_B$ is too low, and traffic-jam-like behaviour are prevented, the plant would also not be able to take up boron, thus, this would be a parameter setting that is biologically irrelevant. Therefore, although our interpretation of the experimental results (Figure 1—figure supplement 2) suggest larger differences between $k_B$ and $k_N$, given that these are indirect estimates that rely on underlying assumptions, this robustness analysis shows that, within a biologically reasonable window, our results hold even when $k_B$ were to be smaller.

Figure 3 with 1 supplement see all

Download asset Open asset

We next wondered if this traffic-jam-like effect dissipates over larger tissues. Our four-cell model is a simplified representation of A. thaliana, which has less cell files between the epidermis and the xylem than most other roots of experimental plant models. However, if one distinguishes the pericycle from the vasculature, a five-cell model would be more appropriate. We show that the variables of the five-cell model present similar dynamics as our default model when BOR transporter regulation is lowered to 400 μM (Figure 4—figure supplement 2A–D). Lowering $k_B$ slightly stabilises the dynamics (as shown in Figure 2—figure supplement 1), which counteracts increasing destabilisation of the constant flow regime when the number of cell files in the tissue increases. Conversely, given that this renders smaller tissues (such as our default four-cell model), more robust against traffic-jam-like behaviour, our default setting should be regarded as a worse-case scenario for such dynamics to occur, rather than being a special case. To further appreciate tissue-size dependencies, and explore the generality of these effects for other plants which might have highly deviating cell file numbers, we gradually increased the in silico tissue from 4 to 10 cell files. Figure 3B–G shows that the phenomenon is conserved irrespective of the number of cells between the medium and the convective xylem flux. Larger transversal tissue segments are able to generate higher wave amplitudes, but in all cases the peak in boron propagates contrary to the nutrient flow with approximately the same speed. The same phenomenon was still observed even when the model contained 100 cell files (Figure 3—figure supplement 1), suggesting that large tissues do not prevent these effects from occurring, but rather increase its likelihood.

In short, when transporter regulation swiftness is sufficiently slow, due to the downregulation of transporters under the control of increased boron concentrations, a traffic jam-like behaviour emerges: a high boron concentration peak appears, correlating with locally lower transport rates, triggering a low-transport high-boron wave that back-propagates from cell to cell in the direction contrary to the transport polarisation direction within the tissue. This occurs independently of the size of the plant tissue under consideration.

We next studied the importance of the relative dynamics of BOR and NIP regulation to trigger these phenomena, that is, which of the processes have to be sufficiently slow.

Under what conditions and parameter values for transporter regulation rates does the traffic jam phenomenon manifest itself? In our previous simulations, we simultaneously decreased the swiftness of both the NIP and the BOR regulation. Oscillations in boron concentration then arose due to a change in stability of both the flow and steady state of boron concentrations. We next probed the system’s stability while varying the regulation swiftness of the NIPs and BORs independently. This was done by analysing the equilibrium of the simplified, A. thaliana inspired, four-cell model. We linearised the ODEs around the equilibrium using a Taylor expansion, and then evaluated the stability of the equilibrium in terms of the largest eigenvalue, for 40,000 different parameter combinations. The phase-portrait that emerges (Figure 4A) shows a large region in which the system becomes unstable and oscillations arise (indicated in red), as well as a region in which sufficient regulation swiftness ensures that the system is stable and oscillations do not develop (indicated in blue). For the stable (blue) parameter settings, any perturbation dampens out (Figure 4B), whereas for parameter settings within the red region, oscillations dominate. The frequency of the oscillations is determined by a relative lack of swift regulation for both transporters, as shown by the colour map in Figure 4. This phase diagram is qualitatively unaffected for the five-cell model at $k_B=$400 μM (Figure 4—figure supplement 2F). Furthermore, we also checked if saturation in the activity of the BOR transporter might obstruct or stabilise the system against the appearance of the oscillations. (The NIPs, being a channel, should be less prone to undergo concentration-dependent saturation effects in permeability.) Contrary to this, we found that such a saturation in permeability of the BOR transporter does not change the likelihood or transporter-response dependency to present traffic-jam-like behaviour, but rather, when they arise due to slow dynamics, the oscillations present much higher amplitudes in concentration variations, with peak boron levels becoming more than twice as high due to the transport saturation (Figure 4—figure supplement 3).

Figure 4 with 3 supplements see all

Download asset Open asset

From this combined analysis it becomes apparent that the parameter values we derived to capture wild type dynamics of transporter regulation (i.e., point 1,1) fall into the stable regime. Hence, the wild type system is robust against perturbations and does not present traffic jam-like behaviour. However, independent small reductions from these values in either NIP or BOR regulation rate can cause the tissue to experience oscillations. Not only does the analysis reveal that both transporters are required to present a rapid response, and that this is independent of the BOR-dependent fluxes becoming saturated, but their combined swiftness is synergistic: if one transporter is extremely fast, the other can be bit slower than what is otherwise needed. However, for any parameter setting, a certain speed in responsiveness needs to present to prevent oscillations.

Traffic jams are commonly experienced by vehicle drivers in urban areas. In terms of human transport, traffic jams are an undesirable outcome as they reduce the throughput through the highway and increase the time to reach the destination. Traffic jams displace backwards in space (contrary to the flow) while giving rise to increased car densities at the specific location of the traffic jam (Kesting and Treiber, 2013). All those properties can also observed in our plant tissue model, namely diminished throughput of boron from the soil to the xylem and backwards propagating pulses of locally increased boron concentrations within cells. Moreover, again analogous to real-life traffic jams, traffic-jam-like behaviour in plant nutrient uptake has undesirable implications in the biological context as well. High intracellular boron is detrimental to plant health as it causes DNA damage and ultimately causes cell death (Sakamoto et al., 2011-09). Plants are therefore likely to have evolved mechanisms to avoid high intracellular boron levels to limit boron-induced damage. Lower throughput across the root implies reduced xylem loading and hence a reduction in the boron absorption efficiency, critical for plant growth at low boron conditions.

To better quantify the magnitude of these effects, we performed a simulation screen to evaluate the expected physiological impact of the transporter swiftness in the form of increased boron levels within the cells due to the instabilities, as well reduction in throughput (see Materials and methods for details). We observe increasing maximum levels in boron concentration for slower regulation swiftness (Figure 5A), as well as higher fold-changes in boron concentration (Figure 5B). This suggests that sufficiently rapid transport regulation is important to reduce the risk of DNA damage when either considering absolute concentrations or the increase over the basal equilibrium boron concentrations. Moreover, under the conditions in which traffic jams and large-amplitude variations occur, we measured a considerable reduction in total throughput through the tissue (Figure 5C). In our default model we assume the BOR transporter-driven flux to be linear with the cellular boron concentration. If we instead consider Michaelis-Menten saturated BOR transport (using a biologically reasonable Michaelis constant — the boron concentration at which the flux becomes half as large due to saturation, as well as the concentration at which the flux reaches its half-maximum value — of $c_B=$1000 μM) the system’s capacity of uptake and throughput only marginally decreases. In stark contrast, the detrimental effects of the instabilities that arise due to slowing down the regulatory dynamics are amplified when BOR saturation is considered (see details of saturation implementation in the caption of Figure 4—figure supplement 3). The impact of the traffic jams becomes much more severe, with peak levels during the concentration fluctuations more than twice as high as the default scenario (compare Figure 2 to Figure 4—figure supplement 3). This is a direct consequence of the cell not being able to effectively efflux boron when intracellular levels spike. Taken all into consideration, we infer that selective pressures operating on the root’s capacity to absorb boron optimally and robustly could result in the system evolving to a regime of rapid regulation of transporters.

Figure 5

Download asset Open asset

Our model thus far consisted of a segment of polarised tissue, of a variable number of cell files, linking the soil to the xylem via a polar flow of nutrients. Consequently, both soil and xylem, although presenting stable characteristics over time, do constitute abrupt boundary conditions: the first sets a constant medium concentration, and the second presents a constant convective flow away from the transversal section. We thus questioned if these boundary conditions are responsible for triggering the traffic jam-like regime that arises under lower rates of transport gene regulation. During the last decades a similar discussion has been taking place regarding what triggers the analogous traffic jams in discrete macroscopic systems such as road networks (Kesting and Treiber, 2013). Although traffic jams are readily triggered by bottlenecks (such as road obstructions, roundabouts, etc.), theoretical models of traffic jam behaviour predicted that the collective systems parameters alone could be sufficient to render the free flow state unstable. It depends on the car density exceeding a critical value as well as how speed relates to the local car density. To prove this theoretical insight, Sugiyama et al., 2008 developed an experiment in which cars were confined to move on a homogeneous circular road. Above a critical car density they observed a transition from a free flow to a traffic jam state, due to the collective effect of the vehicles. They concluded that, although no bottleneck was present, the intrinsic parameters of the system rendered it unstable, such that the smallest of naturally occurring fluctuations was enhanced to disrupt the free flow. This experiment illustrates the non-intuitive notion of traffic jams being generated spontaneously under certain critical parameter regimes.

Inspired by this strategy, and for us to be able to rule out that the origin of the nutrient-flow instabilities within the plant tissue is provided by the soil boron concentration or the xylem flow, we constructed an in silico polarised tissue composed of 6 cells, wrapped up into a ring, that is, the first cell is connected to the last. The ring structure avoids boundary conditions (i.e., bottlenecks) and offers a scenario in which uninterrupted flows (either stable or unstable) can be studied. We allow the dynamics of the polarised tissue to evolve from an initial situation in which concentrations are homogeneous over the whole ring-tissue. When transporter regulation is as swift as in wild type (Figure 6A), a steady flow with no oscillation in the concentrations is observed. When transporter regulation swiftness is decreased to sufficiently low values as to cause unstable flows (as derived from the analysis shown in Figure 4), oscillations arise. After an initial period in which the flows seem constant, random numerical fluctuations in the simulation bring forth a rise of boron concentration in a random cell, which, due to the mechanism of inhibition by boron of the boron transporters, results in local interruption of flows and the back propagation of the boron concentration pulse (Figure 6C,D, Video 1). The frequency of these oscillations is decreased, and the wave broadened, as the regulation speed is decreased (Figure 6C,D). These simulations demonstrate that the unstable flow regime — which we now can conclusively eliminate as being triggered by boundary conditions — continuously self-propagates.

Figure 6

Download asset Open asset

Video 1

Download asset

Our model analysis highlights the requirement for rapid regulation of both BOR1 and NIP5;1 transporters to stabilise boron flux through root cells and minimise risks associated with transient high levels of boron in cells, both considered to be important constraints for the plant. Current experimental evidence shows rapid downregulation of the BOR1 protein (Takano et al., 2005; Takano et al., 2010), and NIP5;1 transcript (Tanaka et al., 2011; Tanaka et al., 2016), although the fine details regarding the molecular mechanisms remain to be elucidated. The modelling however predicts permanently fast regulation, not only when boron levels drop, but also when they remain constant. This is however difficult to prove with techniques that we have used previously.

To experimentally explore the regulatory swiftness of the bidirectional transporter NIP5;1 at the transcript level, we used single molecule RNA FISH (smFISH). In this method ~40 complementary fluorescently labeled oligonucleotide probes are used to visualise RNA at the cellular level (Raj and Tyagi, 2010; Duncan et al., 2017). Initially, we combined probe sets to detect mRNA for both NIP5;1 and the housekeeping gene PP2A (specifically A3 scaffolding subunit of PP2A, At1g13320) and compared abundance and cellular distribution (Figure 7A) (Czechowski et al., 2005; Duncan et al., 2016). Consistent with published smFISH experiments, PP2A mRNA appeared evenly diffused throughout all cells (Figure 7B) (Duncan et al., 2016). In contrast, we observed many nuclear accumulations of NIP5;1 RNA that were accompanied by few cytoplasmic mRNA copies (125 nuclear accumulations in 149 cells) (Figure 7B). Large nuclear smFISH RNA signals have been reported for highly transcribed genes. They are typically accompanied by many copies of cytoplasmic mRNA and considered to represent bursts of transcription, where multiple Pol II associate with a locus (Battich et al., 2015; Duncan et al., 2016).

Figure 7 with 1 supplement see all

Download asset Open asset

Arabidopsis genes that contain introns must undergo pre-mRNA splicing before mature transcripts are translated (Morello and Breviario, 2008-06). Co-transcriptional splicing is common for plant genes and this ensures that introns are removed from the pre-mRNA and rapidly degraded close to the site of transcription (Reddy et al., 2013-10). This system allows intronic smFISH RNA labelling to be used as an effective method to identify loci actively engaged in transcription (Levesque and Raj, 2013; Duncan et al., 2016; Rosa et al., 2016). We used this approach to further investigate NIP5;1 transcription and found 72% of cells with at least one NIP5;1 intron signal. As expected, this was lower than the 84% previously reported for the more highly expressed gene PP2A (Duncan et al., 2016). When we combined NIP5;1 exonic and intronic RNA smFISH probe sets we observed 92.5% of nuclear accumulations co-localised with intron signals ($n=111$) (Figure 7C). This is consistent with the nuclear accumulations representing sites of ongoing RNA production and degradation, rather than separate nuclear storage compartments.

The suggested rapid degradation of NIP5;1 mRNA after transcription was further supported by an independent approach, in which we performed qRT-PCR using probes specific to pre- or mature mRNA to measure mature mRNA/pre-mRNA ratio in root cells. The average mature mRNA/pre-mRNA ratio of NIP5;1 was less than 25% of that of PP2A (Figure 7—figure supplement 1). In accordance with our interpretation of the smFISH results, this again indicates a high degradation rate of NIP5;1 mRNA.

In light of our model predictions, these observations combined point to a highly dynamic sensing system where RNA is turned over at or near the site of transcription, to limit NIP5;1 levels when boron is above a threshold level. We recently demonstrated that NIP5;1 transcript degradation occurs through ribosome stalling, triggering degradation in the cytoplasm (Tanaka et al., 2016). We also identified an upstream sequence that promotes mRNA degradation, but does not trigger ribosome stalling. The data presented in Figure 7 provides evidence of mRNA degradation occurring close to the site of transcription. Together these results suggest that mRNA degradation could occur through two coordinated B-dependent processes; one where cytoplasmic degradation is triggered by ribosome stalling and another localised in the nucleus where RNA is turned over at the site of transcription. Relief of this repression could then ensure rapid boron uptake when required. Our model suggests that such a dynamic mRNA production/decay process could support the necessary swift boron-mediated transporter regulation critical for the evolution of controlled, stable boric acid flows through polarised tissue.

These results also support the notion of a constant turnover (both producing and breakdown the RNA) occurring for regulating NIP5;1 transcript levels. We propose that this ‘wasteful’ process can now be understood in light of constraints of stable flows through polarised tissue. Interestingly, the regulation of the transporters is of a different nature, one predominantly on the level of protein degradation (BOR1), the other occurring on the level of transcript accumulation (NIP5;1).

We conclude that the experimentally observed swiftness of boron transporters’ boron-mediated downregulation can be understood as a necessary mechanism to maintain optimal constant xylem loading rates over time, while also avoiding DNA damage due to oscillations within cells. One could argue that without a substrate-mediated downregulation of the transporters there would not be the need to avoid traffic jam-like behaviour through swift regulatory dynamics in the first place. However, were the system not to present the regulation, the plant would be unable to shield itself from toxic soil boron levels (Figure 8A) and be much less efficient in growing at low soil boron levels (Figure 8B). Thus, some kind of inhibitory regulation necessarily has to be present, to ensure plant plasticity under different environmental conditions. As a consequence, the speed of the transporter downregulation needs to be sufficiently high. We showed that this is an intrinsic requirement, simply because boron transporters are polarly located: The requirement itself does not depend on external boron fluctuations in the medium, number of cell files composing the root tissue, possible saturable activity of the transporter itself nor the magnitude, location or strength of the xylem flow.

Figure 8

Download asset Open asset

Our study has been based on the well-characterised system of directed boron transport, owing to the richness in quantitative experimental measurements regarding BOR and NIP regulation. However, the implications of our derived insights apply to any system in which the dynamics of polar transport of a given substance relies on an inhibitory feedback between the concentration of that substance and its local transport rate. For example, the phosphate and iron transport systems present both polar transporters as well as substrate-dependent regulation, serving as candidates for such phenomena and dynamical constraint avoidance. In the case of the phosphate transport system, transporters localised to the outer side of the epidermis are responsible for uptake (LePT1 and LePT2, Liu et al., 1998), and they have recently been found to carry a domain possibly involved in phosphate sensing which leads to rapid degradation of phosphate transporter in a phosphate-dependent manner (Gu et al., 2016). Similar regulation is also reported for the iron (Fe) transport system, where polar accumulation of IRT1, a major Fe transporter for Fe uptake into roots, is regulated in an Fe-dependent manner (Barberon et al., 2014). Even other transport systems which do not shuttle nutrients might display this sort of instability, such as polar auxin transport, for which stable flows might also require fast response dynamics regarding the activity of auxin importers and efflux facilitators (Robert and Friml, 2009). Indeed, for any polarised tissue in which adaptive regulation of transport levels is necessary we can expect to find intrinsic temporal constraints operating to ensure that steady state flows are maintained and large fluctuations avoided.

In molecular terms, we discovered that although both transporters are predicted to work synergistically and are both required to be effectively similar in regard to rapid regulation, their regulations are biologically encoded in unique ways. We here experimentally focused on the regulatory swiftness of the bidirectional transporter NIP5;1, and found strong evidence that it constitutes a system in which production is maintained high but degradation is rapidly controlled transcriptionally, allowing for the necessarily rapid boron-dependent response. It will be interesting to speculate what underlying reason, if any, has caused the directed exporter, BOR1, which is similarly known and required to be swiftly regulated, to be controlled on a biologically distinct level.

For this work, the analogy with motorway traffic jams was helpful to build an understanding of principles of polarised transport dynamics, albeit both systems obviously deviate in important aspects. Boron concentrations are continuous and can reach arbitrarily high levels, as opposed to (incompressible) cars. Yet, the discrete nature of cells and their regulated polar transporters, combined with the boron-dependent inhibition of the transporters themselves, resembles cars slowing down in response to busy regions along the motorway. Moreover, in both systems highly sensitive responses can prevent jamming (Sugiyama et al., 2008). We therefore propose that the notion of traffic jams can be universally instructive and helpful for biologists considering physiology and regulated growth through polarised transport in plants, as well as for future efforts to enhance plant growth under diverse environmental conditions.

For plant culture, MGRL medium was used (Fujiwara et al., 1992) supplemented with 1% sucrose and solidified with 1.5% gellan gum. For Figure 1, A. thaliana transgenic lines carrying BOR1-GFP (Kasai et al., 2011) and GFP-NIP5;1 #8 (Tanaka et al., 2011) were germinated and pre-incubated for 5 days on normal MGRL medium plate, and then grown for 2 days on MGRL plates with 0.3 μM boric acid. For Figure 1—figure supplement 2, plants were grown for seven days on MGRL medium with the indicated boric acid concentrations. Images were captured with a confocal microscope (FV1000, Olympus, Japan). GFP fluorescence was observed with 473 nm excitation and 510 nm emission. Cell wall was stained with 10 μg/mL propidium iodide aqueous solution for 3 min, and observed with 559 nm excitation and 619 nm emission.

NIP5;1 promoter activity was observed using a transgenic A. thaliana strain carrying NIP5;1 promoter-GUS (strain −2,180$\mathrm{\Delta }$UTR312 #8 in Tanaka et al., 2011). Seedlings were germinated and cultured on MGRL medium plate. Four-day-old seedlings were stained with GUS staining solution: 100 mM Na${}_2$HPO${}_4$ (pH 7.0), 0.1% Triton X-100, 2 mM K${}_4$Fe(CN)${}_6$, 2 mM K${}_3$Fe(CN)${}_6$, 0.5 mg/mL X-Gluc (5-bromo-4-chloro-3-indolyl $\beta$-D-glucuronide cyclohexylammonium Salt, Glycosynth, UK) in a decompressed desiccator for 20 min at room temperature. After rinsed with phosphate buffer (pH 7.0), stained seedlings were mounted in 45°C molten 5% agar (gelling temperature 30–31°C, Nacalai, Japan) aqueous solution and solidified at 4°C for 3 hr. The samples were trimmed into blocks and sliced with a vibrating microtome with thickness of 100 μm. The sliced sections were observed with microscope.

Although we have constructed our model as parsimonious as possible to qualitatively explore what kind of behaviours could emerge from a polarised tissue with transporter regulation, our model requires a few important parameters to be set. We have sought to explore the possible dynamics of our system while staying well within the expected limits of what is experimentally considered plausible in terms of the boron transport system. Below we give a brief description of the data our parameter choices were based on, and, when possible, how reasonable ranges could be established.

The plasmamembrane boron permeability rate ($p$) is set to a maximum rate of 8·10^-2 s${}^{-1}$ and a minimum rate of 4.4·10^-3 s${}^{-1}$, based on the estimation of boric acid membrane permeability by Raven, 1980 and on more recent experimental measurements performed on internodal cells of the charophyte alga Chara corallina (Stangoulis et al., 2001). Note that the plasmamembrane’s permeability rate for boric acid is relatively high compared to other nutrient elements, due to boric acid existing in a non-charged form under physiological pH.

The degradation rates ${\xi }_B$, ${\xi }_N$ represent the effective decay rate of the actual transporter activity, not just of the transporter protein. This includes processes such as its removal from the membrane, its inactivation due to usage, etc. As it amalgamates a range of inactivation/degradation processes, we set it higher than typical protein degradation rates, but lower than typical membrane residence times, using a half life of 9 s. The valid ranges for the production rate of the transporter activities, ${\alpha }_B$ and ${\alpha }_N$, were derived from the degradation rates ${\xi }_B$, ${\xi }_N$, combined with qualitative conditions for the transporter equilibrium levels. It is as yet not possible to obtain these values directly through experiments, as several intermediate steps are involved. Instead we impose that the transporter activity steady state level does not exceed the plasmamembrane’s permeability by a thousand fold:

Combining these conditions we could then extract the valid ranges for ${\alpha }_B$ and ${\alpha }_N$. As explained in the main text, we take ${\displaystyle k_B\simeq 30k_N}$, based on our observations of the sensitivity of BOR and NIP5;1 expression under different boron concentrations (Figure 1—figure supplement 2), further corroborated by earlier data on NIP5;1, as displayed in Figure 1B in Tanaka et al., 2011. We also ran simulations to show robustness of the observed behaviours for smaller variations between $k_B$ and $k_N$ by assuming lower values of $k_B$ (Figure 2—figure supplement 1, Figure 4—figure supplement 2). The parameter $d_B$ determines the maximum fold increase in BOR degradation at very high intracellular boron concentrations. It was set to a 50-fold increase, based on rough assessments of experiments in which plants were transferred to high boron conditions, with the caveat that assessment of the temporal change in functional, membrane-located BOR as a function of the intracellular boron concentration is very challenging. Typical cell sizes (height, $hc$; width, $lc$) and cell wall width ($lw$) were estimated from confocal microscope images of A. thaliana roots. The xylem loading rate $a$ was calculated as follows: Hosy et al., 2003 reported that in mature Arabidopsis the rate of transpiration is 0.1 g water/g FW tissue/hour at 70% humidity under daytime conditions. Given an Arabidopsis FW at 18 days of around 0.2 g/plant, transpiration is 20 mg/hour/plant, or roughly 5·10⁶ μm³/s. Considering a stem xylem cross-sectional surface area being roughly 10⁷ μm², the convective velocity $a$ should therefore be around 5·10⁶ μm³/s / 10⁷ μm² = 0.5 μm/s.

Numerical and analytical calculations were computed using Wolfram Mathematica 10.2. To calculate time development of the model, ODEs were solved numerically using the NDSolve[] function, setting the initial conditions for each variable to zero. The parameters used in the models are shown in Table 1. The solutions of the ODEs were visualised as a line plot or kymograph, using the Plot[] or ListDensityPlot[] functions.

To evaluate the physiological impact (maximum concentration and amplitude of cytosolic boron, and throughput towards shoots) of transporter regulation swiftness, numerical simulations were repeated for combinations of ${ϵ}_B$ and ${ϵ}_N$, which were independently varied between 0.01 and 2, at intervals of 0.01. For each simulation, minimum and maximum boron concentration for each cell and average throughput through the cell row between 24 to 48 h simulation time were recorded using the ‘EventLocator’ method of the NDSolve[] function.

For the ring model, the boundary conditions in the original model were removed by connecting the last cell in the tissue-segment to the first one. This was done by using a tissue composed of six cells and then introducing a cell wall between cell 6 and cell 1. The boron concentration in the novel cell wall is referred to as $W_{6-1}$. This generates a six-cell ring model. Because within a ring-model the total amount of boron is conserved, introducing an additional ODE for this cell wall component effectively means an overdetermination of variables, causing numerical issues. Hence the value of the new variable was instead defined as follows:

where $B_{total}$ is the total amount of boron in the whole system, which is constant over time as there is no input into nor outflow from the system. The initial boron concentration of all cells and cell walls in the simulation were set to 300 μM, and $B_{total}$ was determined accordingly.

The local stability analysis was performed as follows: First, all possible equilibria for the set of ODEs were found using the NSolve[] function. Next, analytical descriptions of the Jacobian matrix of the ODEs around the equilibria were derived using the D[] function. Finally, the stability of the system for varying values of ${ϵ}_B$ and ${ϵ}_N$ was determined according to the signs of the eigenvalues of the Jacobian matrix. (Please note that changing ${ϵ}_B$ and/or ${ϵ}_N$ does not change the number of equilibria or their value, only their stability can change.) The system is stable when the real parts of all the eigenvalues are negative, and unstable in any other case. The imaginary part of the largest eigenvalues was then used to estimate the oscillatory period in the unstable region, based on the following conversion:

To compare the predicted oscillatory period close to the unstable equilibrium with the actual oscillatory period of the system, the oscillatory period was also calculated numerically (Figure 4—figure supplement 1). For certain parameters individual variables can present rather complicated patterns, for example temporarily maximum values, followed by a slight decrease and then further rise before decreasing again. We therefore used the following algorithm to determine the oscillatory period: The ODEs were numerically solved using the NDSolve[] function. The ‘EventLocator’ method was then used to collect all extrema for three variables, $C_2$, $C_3$ and $N_2$, for the period between 24 and 72 hr. Using the first set of extrema as a reference point, the second and potentially later sets of extrema were then compared against the first set. The first case for which all differences were within ±0.1% for all three variables was then labeled as being identical to the reference point, and the corresponding time interval reported as the oscillatory period.

The probes were designed using the online program Stellaris™ Probe Designer version 2.0 and ordered from LGC Biosearch Technologies (Petaluma, CA). For probe sequences see Table 2.

Col-0 seeds were stratified for 2 d at 4°C and then germinated and grown under 16/8 hr light conditions at 20°C on MGRL plates (30μM boric acid). smFISH was carried out on seedlings as described by Duncan et al., 2017. Briefly, 5d old seedlings were fixed for 30 min in 4% paraformaldehyde. Root squash samples were prepared on slides. Tissue permeabilisation was achieved by immersing the samples in 70% ethanol for a minimum of one hour. Two washes were carried out with wash buffer (10% formamide and 2x SSC). 100 mL of hybridisation solution (0.3 M Sodium chloride and 30 mM tri-Sodium citrate dihydrate, referred to as ’2x SSC’, 10% dextran sulfate and 10% formamide) containing probes (final concentration 25 nM) was then left to hybridise at 37°C in the dark overnight. Each sample was washed twice with wash buffer following hybridisation buffer removal with the second wash left to incubate for 30 min at 37°C. Each slide was then incubated with 100mL 4’ 6-diamidino-2-phenylindole dihydrochloride (DAPI, Sigma, St. Luis, MO) (100 ng/mL) at 37°C for 30 min. The DAPI solution was removed before 100 μL 2x SSC was added and then removed. Samples were then incubated for 2 min at room temperature with 100 μL GLOX buffer (0.4% glucose in 10 mM Tris, 2x SSC) before being mounted in 100 μLof GLOX buffer containing 1 μL of glucose oxidase (#G0543, Sigma) and 1 μL catalase (Sigma) under 22 mm x 22 mm No. 1 coverslips (Slaughter, Uppminster, UK) and sealed by nail varnish. A Zeiss Elyra PS1 inverted microscope was used for imaging. A x100 oil-immersion objective (1.46 NA) and cooled EM-CCD Andor iXon 897 camera (512$\times$512 QE>90%) were used to obtain all images. (Quasar^®570 probes: excitation by 561 nm laser, signals detected between 570–640 nm. Quasar^®670 probes: excitation by 642 nm laser, signals detected between 655–710 nm. DAPI: excitation by 405 nm laser, signals detected between 420–480 nm.) For all experiments a series of optical sections were set up with z-steps of 0.2 μm. All smFISH images were de convolved with Huygens Professional version 16.05 (Scientific Volume Imaging, the Netherlands, http://svi.nl) using the QMLE algorithm with Good signal to noise settings and 50 iterations. Maximum Z projections in Figure 7 were performed using Fiji (an implementation of ImageJ, a public domain program by W. Rasband available from http://rsb.info.nih.gov/ij/).

Col-0 seedlings were germinated and grown on MGRL plates containing 30 μM boric acid for eight days. Total RNA was extracted from whole roots of about 25 seedlings with four replicates, using NucleoSpin RNA Plant (MACHEREY-NAGEL). To reduce contamination of genome DNA, 500 ng of prepared RNA was subjected to DNase treatment with RQ1 RNase-Free DNase (Promega), following the manufacture’s protocol. Reverse transcription was performed with Super Script III reverse transcriptase (Thermo Fisher Scientific) with random primers. To compare the abundance of cDNAs with different sequences, absolute quantification by realtime PCR was conducted, using PCR-synthesised DNA fragments with known concentration as standard templates. To detect mature or pre-mRNA specifically, exon-exon or intron-exon probes were designed (Figure 7—figure supplement 1 and Table 3). Realtime PCR was performed with Thermal Cycler Dice® Real Time System and SYBR Premix Ex Taq II (Tli RNaseH Plus) (Takara), following the supplied standard shuttle 2-step PCR protocol. As standard DNA templates, DNA fragments corresponding to cDNA of mature- or pre-mRNA of NIP5;1 and PP2A were amplified by PCR with Prime Star GXL (Takara), using Col-0 cDNA or genome DNA as a template. The primers used are shown in Table 3. The PCR products were subjected to electrophoresis and the target fragment was purified from the gel slices with QIAquick Gel Extraction Kit (Qiagen). The concentration of the purified DNA was estimated by NanoDrop 1000 (Thermo Fisher Scientific). Tenth dilution series of 1 pg/μL~1 ag/μL template DNA were applied for standard curve. The copy number of each standard DNA solution was normalized by realtime PCR $C_t$ values of the exon-exon proves which amplify both mature- and pre-mRNA sequences (Figure 7—figure supplement 1 and Table 3). As a negative control, in parallel, after the total RNA extraction, the same samples were processed with the same procedure except reverse transcription, confirming that realtime PCR signals from contaminated genome DNA does not affect the final results.

Video 2

Download asset

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

This work has been supported by the UK Biological and Biotechnology Research Council (BBSRC) via grant BB/J004553/1 to the John Innes Centre and by the Japan Society for the Promotion of Science via Grants-in-Aid for Scientific Research (No. 15J11021 to NS and No. 25221202 to TF). Experiments were partly financed by the OpenPlant BBSRC/EPSCR Synthetic Biology Research Centre (Grant BB/L014130/1).

© 2017, Sotta et al.

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