Regulatory feedback response mechanisms to phosphate starvation in rice

Phosphorus is a growth-limiting nutrient for plants. The growing scarcity of phosphate stocks threatens global food security. Phosphate-uptake regulation is so complex and incompletely known that attempts to improve phosphorus use efficiency have had extremely limited success. This study improves our understanding of the molecular mechanisms underlying phosphate uptake by investigating the transcriptional dynamics of two regulators: the Ubiquitin ligase PHO2 and the long non-coding RNA IPS1. Temporal measurements of RNA levels have been integrated into mechanistic mathematical models using advanced statistical techniques. Models based solely on current knowledge could not adequately explain the temporal expression profiles. Further modeling and bioinformatics analysis have led to the prediction of three regulatory features: the PHO2 protein mediates the degradation of its own transcriptional activator to maintain constant PHO2 mRNA levels; the binding affinity of the transcriptional activator of PHO2 is impaired by a phosphate-sensitive transcriptional repressor/inhibitor; and the extremely high levels of IPS1 and its rapid disappearance upon Pi re-supply are best explained by Pi-sensitive RNA protection. This work offers both new opportunities for plant phosphate research that will be essential for informing the development of phosphate efficient crop varieties, and a foundation for the development of models integrating phosphate with other stress responses.

compartmental phosphate levels -Cytosolic Pi (CytoPi) (in seconds and minutes) and the induced genetic/molecular response (hours and days) occur at different time-scales.
In the model, both low-affinity and high-affinity phosphate transporters are assumed to be embedded in the cell membrane, corresponding to all families of Pi transporters. The constitutively expressed low-affinity transporters account for some basal Pi uptake under all conditions. The rate of phosphate utilization is assumed to be solely dependent upon availability of Pi in cytosol. The external Pi concentration in the model is considered to remain constant under typical and deficient conditions at 200 μM and 0.0001 μM respectively. Given the upper limit for the concentration of cytoplasmic phosphate (i.e. 25mM) [7] and the proportion of it stored in the vacuoles (~80%) [8,9], the CytoPi level is set at a steady state value of 5000 μM (i.e. 5mM) under typical hydroponic external Pi conditions, i.e. 200 μM.
The abundance of each mRNA species is used as a proxy for its protein concentration, assuming that mRNA concentration is at quasi steady state. Because the molecular regulation of Pi uptake is largely conserved between rice and Arabidopsis [1], the halflives for all the transcripts in the model correspond to those of Arabidopsis [10]. Published and qRT-PCR results shows that PHR2 mRNA levels (and hence total protein) remain roughly constant, in both typical and low Pi conditions. PHR2 is defined as the level of unbound PHR2 available for sumoylation by SIZ1. SUMOylation of PHR2 is countered by de-SUMOylating enzymes, potentially SUMO proteases [11,12], whose activity is assumed to be constant.
The rate of PHO2 mRNA synthesis is also assumed to be constant because there is no evidence to the contrary. The binding of IPS1-mR399 is weaker than PHO2-miR399 binding because of the 3bp insertion in the middle of miR399 binding site in IPS1, while the dissociation rate of the IPS1-miR399 (IMC) complex is assumed to be lower than their association rate.

Equations
The model is a system of ten coupled non-linear ordinary differential equations (ODEs).
The corresponding model variables are: the concentration of active SIZ1 protein; the concentration of free un-sumoylated PHR2 transcription factor; the concentration of sumoylated PHR2 transcription factor; the concentration of miR399 RNA; the concentration of PHO2 mRNA; the concentration of PHO1 protein; the concentration of membrane-bound high-affinity PHT protein; the concentration of IPS1 RNA; the concentration of IPS1-miR399 complex (IMC); and the concentration of CytoPi in roots.
The rate of endogenous production and degradation are respectively denoted by "m" and "d", along with relevant subscripts and square brackets denote concentrations. The equations are presented below.
The active SIZ1 protein is repressed by cytosolic phosphate (CytoPi) with the inhibition constant k1 and Hill co-efficient n1 (Equation 1). Under typical Pi conditions, the transcription factor PHR2 remains bound to Pi-stabilised SPX4 protein [13]. For simplicity, PHR2 is assumed to interact directly with CytoPi (Equation 2). The free unbound PHR2 is sumoylated with Michealis constant k3 and eventually, de-sumoylation with rate constant k4 (Equation 3). The bulk translocation of miR399 from shoot-to-root [ 1] =

Known, assumed and calculated parameters
Values of 9 parameters, mainly the degradation rates, were obtained directly from the literature [10]. The binding constant for miR399 and PHO2 (k6) is assumed to equal 0.009, on the basis of the work published by [14]. The Hill co-efficient (n) was set to 2 in keeping with many models of gene regulation while Hill co-efficient for IPS1 (r) was set to 4 to achieve the observed delayed and steep increase in IPS1 levels. The inhibition coefficient for SIZ1 (k1), PHO1 and PHT were set to unity, in order to optimize the number of parameters for estimation. The value of m5, k13 and U were calculated as following, using the steady state levels of the respective variables.
Calculation for the production rate of PHO2 (m5) data points [6]. As fold change datasets are used, no units were assigned to the parameters, except for the time in hours (h). The initial values of the variable were set to zero, with the exception of PHO2 and CytoPi, which were initially set to 1 and 5000, respectively.
Model fitting was conducted using a software named MONOLIX that incorporates the suite of parameter inference techniques, detailed in the following section.

Theory
The qRT-PCR data correspond to overall values for a population of multiple root cell

Algorithm
Initially, the individual parameters were simulated using an MCMC approach in the expectation step of SAEM algorithm. These individual parameters were then used to compute a stochastic approximation of the conditional expectation of the log-likelihood of the complete data. Subsequently, the complete log-likelihood was maximized to obtain the updated estimates of the population parameters. In practice, SAEM generates 1 to 5 random samples per individual per iteration to allow for an efficient and rapid convergence toward the solution. This "burn-in" period typically requires 200 iterations.
Thereafter, the program accumulates random sample results for the next iterations, to estimate the population means and inter-individual variances.

Model adjustment
Given that IPS1 RNA is partially complementary to miR399 (with 3 base mismatch), the binding affinity of miR399 toward IPS1 would be relatively weaker than that to PHO2.
However, the release of IPS1 from the IPS1-miR399 complex would be slower than its binding. Considering such a relationship, values of the parameter k7 (the binding constant for IPS1 and miR399) and k8 (Dissociation constant for IPS1-miR399 complex) were indirectly deduced as following: The above algebraic equation were incorporated in the model fitting to first estimate Y1 and Y2, which were then further used to calculate the value for k7 and k8. These parameters (i.e. Y1 and Y2) were purely included to aid the estimation and are related to the known parameter k6 (Binding constant for miR399 and PHO2).

MONOLIX settings
The structural model (i.e. the Pi-regulation model) encoded in MLXTRAN format in a text file (Laveille, 2014). This model and the datasets were loaded onto the MONOLIX platform. Screen shot of the MONOLIX GUI platform is presented in SI.1. Figure 13. The GUI platform is divided into four frames: the data and model, the initialization, the algorithms and the results. Each of these frames requires some input for the precise estimation of parameters. In this estimation exercise following settings were used. In the data and model frame, the 'Combine1' observation model was selected for the model variable to which data are fitted. In the initialization frame, the initial value and standard deviation for the random effect for each parameter were entered. In the algorithm frame, the default Simulated Annealing option was unselected, the number of iteration (K1 and K2) and chains were respectively set to 800, 200 and 50. A new seed was generated for every run, and for other options the default settings were used. In the result frame, linearization was selected for standard errors and log-likelihood, while the conditional mode was selected for individual parameter options. Finally, estimation can be started using the run button given in the top center. Details underlying all these options are well defined and explained in the tutorial available on the MONOLIX webpage (lixoft.eu/products/monolix/documentation).

Estimation settings
Due to the sparsity of dataset for parameters to be estimated, the inference was done by multiple sequential runs as follows. In the first run, the initial values of the parameters and their corresponding standard deviations were set to unity, and the program was executed.
The resulting estimates were than assigned as the initial values for the second run. This cycle was repeated until the best fits were achieved. The algorithm initially struggled to find the range for certain parameters that would result in an acceptable fit to the extremely high IPS1 levels. Using the observation while manually exploring the parameter space, the initial values for m8 and k10 were manually set to 1000 and 70 with standard deviation of unity. These adjustments gave appropriate estimates with acceptable fits to the observed data. The list for the initial values and estimates for the parameters and corresponding standard deviation from each run is given in Supplementary. Table 4.

Estimated parameters for PiOM
The estimated mean values of the all parameters their corresponding standard deviation are listed in Supplementary. small. This is due to lack of data for their respective variables. Perhaps, for a similar reason, the values for k2 and k4 were poorly estimated with higher standard deviation than their mean, but due to their small magnitudes they do not affect the model simulations.

Hypotheses models to explain early PHO2 dynamics in response to Pi-stress
Given the poor prediction of PHO2 profile by PiOM model, five hypotheses are considered as potential mechanistic explanations for the observed PHO2 dynamics in response to Pi-starvation. The first four hypotheses assume some unknown regulator Z acting as (i) a Pi-dependent transcriptional activator of PHO2 (PdTA) -eqs. 11

IPS1 protection hypothesis model
Two strike features concerning the dynamics of IPS1 have been observed. The first is its extreme and unremitting elevation following 21 days of Pi stress and the second is its sudden drop within 24 hours of Pi repletion. The exact reason underlying these features of IPS1 dynamics is vaguely understood. Though the model simulation is somewhat able to explain IPS1 and PHO2 dynamics in response to Pi-starvation. However, the simulation incorrectly predicts the observed drastic drop in the IPS1 and PHO2 level in response to Pi re-supply. This suggests that there must be some extra level of regulation concerning IPS1 and PHO2, which is not correctly represented in the current models.
Compared to other, RNA protection hypothesis was identified as the most plausible explanation for the observed elevation in the level of IPS1 in response to Pi-stress and its sudden drop upon Pi re-supply. Using the PiOM and PsTR models, the IPS1 protection hypothesis was tested by modulating the degradation rate of IPS1 and IMC as the function of CytoPi and examining Pi repletion conditions after a period of phosphate stress.
Furthermore, the magnitude of the IPS1 degradation rate in the RNA Protection (RP) version of the models was altered ( 7 ′ ) and is defined as the quotient of original degradation rate of IPS1 (d7) and the steady-state initial CytoPi concentration under normal Pi condition (CytoPi, at time = 0), see equations (18 -21). In these new models, the degradation rates of IPS1 and IMC will decrease in response to low Pi, allowing IPS1 to accumulate and its reversible interaction with miR399 to form more IMC complexes.
However, the degradation rates of IPS1 and IMC rapidly return back to normal upon Piresupply. The sudden degradation of IMC will thereby cause a rapid increase in the pool of miR399 and consequent short-term decline in PHO2.
PiOM-RP and PsTR-RP model were fitted to the Pi-stress data set and parameters were reestimated. In this round of parameter estimation, previously estimated parameters in Supplementary Tables 3 (for PiOM) Figure SI3.1.
The integrated response (i.e. total sensitivity of the model output across the time-course with respect to the variation of parameters) showed, as expected, that most variables were mainly sensitive to their synthesis and degradation rates. However, most variables (most especially SUMOylated-PHR2) are sensitive to the parameter (m2) -the rate of PHR2 production, Figure S1. As PHR2 is central to the regulation of the phosphate-starvation response, its rate of synthesis is bound to affect various downstream variables. The concentration of CytoPi is most sensitive to the rate of its internal utilisation (U), which had no effect on any of the other variables.
Robustness is the property of a system that maintains its function(s) under perturbations.
Robustness analysis was performed to investigate the robustness of the model variables against the total parameter variation. The Robustness-coefficient identified from this analysis was recorded for all the variables from respective output plots and were replotted together as the bar-chart, Figure S2. For more details about this analysis can be found in [17]. This analysis identifies CytoPi as the most robust variable underlining the fact that the system always tries to keep its concentration constant. The high-affinity phosphate transporters (PHTs) are the next most robust variable, reflecting system's need to use these transporters only when necessary. On the other hand, IPS1 is the least robust relatively more robust upon the inclusion of RNA protection in the models. However, there is no or negligible effect of RNA protection on the robustness of other variables. It also worth noting that IPS1 shows a large variation in its level between different rice varieties upon Pi stress [18].