Modeling reveals posttranscriptional regulation of GA metabolism enzymes in response to drought and cold

To understand how GA distribution is regulated, we constructed a cell-based model that describes GA biosynthesis, degradation, and dilution (due to cell growth) within the maize leaf growth zone (Fig. 1 A and B). The model exploits the simple linear leaf geometry and represents the leaf as a single file of cells. Based on a stable leaf elongation rate during the first 5 d after emergence (31), we considered leaf growth to be in steady state. The model integrates cell growth and division rates from experimental measurements (32) (Fig. 1 C–F). Cell length decreases slightly close to the base of the leaf before increasing with distance along the growth zone (Fig. 1C). Velocity is zero at the base of the leaf and increases with distance along the growth zone (Fig. 1D). Using these data, we calculated the relative elongation rates (Fig. 1E), which follow a roughly bell-shaped curve, and the cell division rates (Fig. 1F), which show a bell-shaped curve spanning the DZ. The model also integrates the increase in the leaf’s cross-sectional area along the growth zone to accommodate increases in cell volume when simulating dilution (33). Consistent with the approximate doubling of both the width and thickness of the leaf (33), volumetric quantifications showed that the leaf cross-section increases more than fivefold across the growth zone (Fig. 1G). With these growth dynamics (Fig. 1 C–F), the maize leaf growth zone is represented by a file of approximately 1,400 cells, with ∼700 cells in the DZ and ∼700 cells in the EZ. These growth dynamics are used to simulate dividing and growing cells (SI Appendix).

In the model, we incorporated the subcellular structure of the cells. Within the DZ, cells predominantly contain a nucleus and cytoplasm; we assumed that the nuclear volume is constant and equal to 50% of the cell volume at the most basal position (noting that the nuclear volume stays constant in the virtual absence of endoreduplication in maize leaves [31]) and thus that cell growth in the DZ occurs due to cytoplasmic expansion. Within the EZ, growth occurs primarily by rapidly increasing vacuolar volume. Based on cell length and cross-sectional area distributions (Fig. 1 C and G), we calculated that cell volume increases ∼12-fold over the EZ, which, under the assumption of vacuolar expansion (with no increase in the nucleus or cytoplasm volume), results in the volume of the cells' vacuole being approximately 92% of the cells’ volume when they enter the mature zone. This value agrees with previous suggestions that the vacuole takes up 90 to 95% of the cell volume (34).

GA biosynthesis involves a series of oxidation steps, converting the precursor, geranylgeranyldiphosphate, to the bioactive GA₄ and GA₁ (35, 36). GA biosynthesis has been shown to be predominantly regulated at the later steps of this pathway (37), whereby GA₅₃ and GA₁₂ are converted to bioactive GAs (35, 36). Focusing on the pathway that leads to the more prevalent bioactive GA in maize, GA₁ (9) (Fig. 1B), we simulated GA biosynthesis and degradation within each cell: GA₅₃ undergoes a series of oxidation steps mediated by GA20ox to produce GA₂₀, which is converted to the bioactive GA₁ by GA 3-oxidase (GA3ox) (35). GA2ox degrade the bioactive GA₁ and precursor GA₂₀ to GA₈ and GA₂₉, respectively, which are in turn converted by the GA2ox to their catabolite forms (35). We represented these reactions by a system of ordinary differential equations (ODEs) for the metabolite, enzyme, and complex concentrations: Each step was modeled using the law of mass action by assuming the GA metabolite first binding to the enzyme with a reversible reaction, and the resulting complex then dissociating into the next GA metabolite in the pathway and the enzyme (38). We assumed that enzymes are translated at a rate proportional to the transcript level.

The reactions involved in the GA metabolism pathway downstream of GA₅₃ occur in the cytoplasm (35, 36), and we assumed that the enzymes and complexes are only present in this compartment. Data in Arabidopsis suggest that GA metabolites are also present in the nucleus and the vacuole (39). In absence of analogous data in maize, we assumed this to be similar in maize, hence within the model assuming equal metabolite concentrations throughout the cell.

Prescribing the growth dynamics and distributions of the GA₅₃ concentration and GA20ox, GA3ox, and GA2ox transcript levels, the cell-based model could be simulated to predict the distributions of the downstream metabolites, enzymes, and complexes. The spatial distributions of the GA₅₃ concentration and GA20ox, GA3ox, and GA2ox transcript levels are upstream inputs, each of which were represented by a sum of b-spline functions (40) with coefficients that were estimated using the experimental data as part of the model fitting (SI Appendix).

To summarize, we developed a cell-based model that describes GA metabolism and dilution within the maize leaf growth zone; the key assumptions behind this cell-based model (described above) are compiled in SI Appendix, Table S1.

The cell-based model comprises 17 ODEs for each cell in the growth zone, which for 1,400 cells results in a system of 25,200 ODEs, which can be simulated until they reach a steady state. For a given parameter set, the cell-based model took several hours to run to predict the steady-state distributions, making detailed parameter surveys impractical. To estimate model parameter values that enable the model to reproduce the experimental data, we needed to derive a reduced model to reduce the simulation time (Fig. 2A). We derived a continuum description of the cell-based model, considering quantities in terms of distance from the leaf base. We further reduced the model by assuming that the ratio between the enzyme concentrations and metabolite concentrations are small (an approximation typically taken when modeling enzyme reactions [41] and shown previously to be appropriate for GA20ox-mediated oxidation [42]). The resulting reduced model involved a system of six ODEs in terms of distance from the leaf base for the concentrations of GA₄₄, GA₁₉, GA₂₀, GA₁, GA₂₉, and GA₈ (highlighted with green stars in Fig. 1B) that depend on eight oxidation rate constants (one associated with each oxidation step, which encompasses the translation rate, binding rates, and enzyme activity), the four input functions (representing the spatial distributions of GA₅₃ concentration and GA20ox, GA3ox, and GA2ox transcript levels; pink stars in Fig. 1B), and the prescribed growth dynamics.

To aid clarity, we provide a summary of the model assumptions underlying the reduced model in SI Appendix, Table S1. As described in the text below, using the reduced model, we were able to estimate the reduced model parameters for a given experimental dataset, and therefore all model results presented in Figs. 2–5 were created by simulating the reduced model.

To test whether the reduced model could represent our observations, we initially parameterized the reduced model using published experimental measurements of metabolite and transcript levels within 12 leaf segments along the maize leaf growth zone, fitting the reduced model parameters independently to data from B104 (9) and B73 inbred lines (30). Prior to fitting, we converted the metabolite measurements (in ng/gDW) (SI Appendix, Figs. S4 and S5) to nM concentrations (SI Appendix, Figs. S6 and S7). After the conversion, the spatial metabolite distribution profiles were globally similar to the original data but could differ in detail; for example, in the nM concentration profile, the peak GA₁ was slightly closer to the leaf base than in the corresponding GA₁ measurements (SI Appendix, Fig. S8). Based on these data, we estimated parameters by minimizing a weighted sum-of-squares criterion (detailed in SI Appendix, section 2.3.4). With the estimated parameters, the reduced model showed a reasonable agreement with the experimental measurements (Fig. 2 B–K and SI Appendix, Fig. S10) and faithfully reproduced the peak in the cytoplasmic GA₁ level within the DZ (Fig. 2I).

Using the reduced model enabled us to estimate the model parameters (i.e., the rate constants in the metabolism network); thus, we were able for the first time to assess the relative impact of individual cellular and subcellular processes on the established GA₁ distribution. Removing either the presence of dilution (Fig. 3A) or the presence of cell movement (Fig. 3B), or the presence of both dilution and cell movement (Fig. 3C), had little effect on the predicted GA distributions. Similar results were obtained for B104 (SI Appendix, Fig. S10). The influence of dilution and cell movement on the GA distribution depends on the magnitudes of the rate constants—with the estimated rate constants, the metabolism network quickly reaches an equilibrium within each cell so that dilution and cell movement are slower processes that have little effect on the predicted GA concentrations (SI Appendix, Fig. S11, which shows how dilution and cell movement have an effect on the GA₁ distribution if the rate constants are smaller). We conclude that dilution and cell movement have only minor effects on the GA distributions and that the GA₁ distribution is predominantly determined by the spatial variations in metabolism.

The estimated parameters provide insights into the mechanisms that determine the distributions of the GA metabolites and differences between B73 and B104. The parameter estimates obtained (SI Appendix, Table S2) suggest that for B73 the GA2ox-mediated degradation rate of the precursor GA₂₀ is small but that there is a faster GA3ox-mediated conversion of GA₂₀ to the bioactive GA₁, which explains the low GA₂₀ concentrations observed. In contrast, for B104, the degradation of GA₂₀ is fast, whereas the conversion of GA₂₀ to GA₁ is slower. To test this model prediction experimentally, we therefore directly compared the rates of GA2ox and GA3ox in the growth zone of B73 and B104 leaves. In agreement with the model prediction, these data revealed that in B73, GA3ox enzyme activity (producing GA₁) is consistently higher than GA2ox activity (producing GA₂₉; Fig. 3D). Moreover, as predicted by the model, in B104 the inverse situation occurs (Fig. 3E).

As one may expect, doubling the GA20ox-mediated oxidation rates increased the predicted GA₁ concentrations (Fig. 3F and SI Appendix, Fig. S12), whereas doubling the GA2ox rates decreased the predicted GA₁ concentrations (Fig. 3F and SI Appendix, Fig. S13), although the qualitative features of the GA₁ distribution remained the same. Varying the oxidation rate associated with the GA3ox-mediated step had little effect on the GA₁ predictions (Fig. 3F and SI Appendix, Fig. S14; doubling the GA3ox oxidation rate increased the rate at which GA₂₀ is converted to GA₁, but this also decreased the GA₂₀ concentrations: at quasi-steady state these processes canceled each other out, resulting in little effect on the GA₁ distribution).

The reduced model also enabled us to investigate the importance of the spatial distributions of the GA₅₃ metabolite and enzyme transcript levels. With constant GA₅₃, the predicted GA₁ formed only a small peak in the DZ and increased as cells left the growth zone (Fig. 3G and SI Appendix, Fig. S15), whereas with constant enzyme levels, the GA₁ peak in the DZ was less pronounced (Fig. 3H and SI Appendix, Fig. S16). We conclude that the spatial variations in both GA₅₃ and enzyme transcription are essential to create the GA₁ distribution that underpins growth regulation.

We next set out to test whether the reduced model could explain the effect of experimental perturbations on GA metabolism, distribution, and leaf growth. We first studied the effects of overexpressing the AtGA20-oxidase1 biosynthesis enzyme (UBI::GA20-OX-1), which enhances bioactive GA levels and growth in both Arabidopsis (43, 44) and maize by increasing DZ size (9) (Fig. 4A). To investigate how overexpressing AtGA20-oxidase1 affects the metabolism dynamics, we simulated the UBI::GA20-OX-1 dynamics by including an additional enzyme, AtGA20ox, in the reduced model and incorporating terms representing the rate at which the AtGA20ox enzyme mediates the three oxidation steps: GA₅₃ to GA₄₄, GA₄₄ to GA₁₉, and GA₁₉ to GA₂₀. Initially, we assumed that for each of these steps the native and heterologous GA20ox enzyme transcripts mediate the same rate of metabolite oxidation, and we tried to fit the reduced model to metabolite and transcript data from wild-type and UBI::GA20-OX-1 (9). With this assumption, we found that the reduced model could not recapitulate the spatial distributions of transcript and metabolite levels observed experimentally (SI Appendix, Fig. S17); the reduced model suggests that for the downstream GAs to be higher in UBI::GA20-OX-1 requires GA₅₃ to also be higher, which is not reflected in the experimental data (Fig. 4 B and G–L).

We solved this conundrum by allowing each of the three GA20ox-mediated oxidation rates to be different between the native and heterologous enzymes, which led to reasonable agreement between the reduced model and data (Fig. 4 B–L). Considering the estimated parameters (SI Appendix, Table S4), the estimated conversion rate of GA₅₃ to GA₄₄ is ∼20 times higher for AtGA20-oxidase1 and the conversion rates of GA₄₄ to GA₁₉ and GA₁₉ to GA₂₀ are approximately double that of the native enzyme. This explains why downstream GA concentrations were higher in the overexpression line (Fig. 4 G–L) despite GA₅₃ concentrations being lower (Fig. 4B). These differences are likely due to differences in translation efficiency, protein degradation, or enzyme activity between the native and heterologous enzymes. The reduced model shows that the differences in GA20ox activity result in the GA₁ concentration having a higher maximum but decreasing to a similar level at the boundary between the DZ and EZ (at approximately 18 mm and 25 mm from the leaf base for wild-type and UBI::GA20-OX-1, respectively) (Fig. 4J), consistent with the GA₁ distribution controlling the DZ size via a threshold mechanism (i.e., the transition to the EZ occurring where GA₁ levels decrease below a threshold value) (9). We conclude that the reduced model recapitulates published data and identifies details in the molecular regulation of GA metabolism, such as differential specificity in the activities of the native and heterologous gene product.

Next, we used the reduced model to determine whether and how GA distributions are affected by environmental stress and if this could explain the growth response. We applied the reduced model to published experimental data involving drought conditions (30) and newly collected data in cold conditions (31). Both stresses reduced the leaf elongation rate by 20 to 30% (30, 31) (Fig. 5 A and B), but this was due to different underlying cellular behaviors: a reduction in DZ size in drought conditions (30, 45) (Fig. 5C and SI Appendix, Table S3) and a reduction in division and elongation rates in cold conditions (31) (Fig. 5D and SI Appendix, Fig. S2).

Measured metabolite levels (in ng/gDW) showed GA₁ levels are reduced in both drought and cold conditions (Fig. 5E). We first converted these measurements to concentrations (nM), taking into account the cross-sectional area. In drought, cross-sectional areas are similar to those of the control (Fig. 5F), so that as for the GA₁ measurements, GA₁ concentrations (in nM) are also reduced (Fig. 5G). However, in cold, cross-sectional areas are much lower (Fig. 5F), resulting in GA₁ concentrations that are in fact similar to those in control conditions (Fig. 5G). We concluded that drought conditions, but not cold conditions, reduce the GA₁ concentrations. Measured GA₅₃ levels (in ng/gDW) were little affected by drought or cold (SI Appendix, Fig. S4), resulting in GA₅₃ concentrations that were increased in cold (once converted to nM), suggesting that cold affects the pathway upstream of GA₅₃ (Fig. 5H). These observations illustrate the importance of conversion to nM concentrations when interpreting metabolite measurements.

Measurements of transcript levels revealed that drought reduces GA20ox levels and increases GA2ox levels (Fig. 5 I and K), whereas cold reduces GA2ox levels (Fig. 5K). We first tested whether these perturbed enzyme transcript levels cause the observed metabolite distributions, assuming identical oxidation rate constants, by fitting the reduced model to the control, cold and drought transcript and metabolite data. With this assumption, we were unable to reproduce the metabolite distributions (SI Appendix, Fig. S18). The predicted difference between GA₁ in control and drought conditions was much less than observed, whereas GA₁ was predicted to be higher in cold than in control conditions, again in contrast to the data.

To resolve this discrepancy, we hypothesized that drought and cold regulate the GA pathway by additional mechanisms (e.g., translation, protein stability, or enzyme activity). To test this theory, we fitted the reduced model to the data, allowing the oxidation rate constants for either the GA20ox-mediated steps, the GA3ox-mediated step, or the GA2ox-mediated steps to be different in cold and drought. To select among the resulting 16 possible cases (SI Appendix, Table S5), we used the Akaike Information Criterion (AICc) (46), a statistical measure that assesses the goodness of fit while penalizing model complexity by taking into account the number of model parameters. Fitting the model and calculating the AICc in each case provided a means to select among the possible cases. Considering the AICc values obtained (SI Appendix, Table S5), our results suggest that the rate constants governing the GA20ox-mediated steps are perturbed in drought and that the rate constants governing the GA2ox-mediated steps are perturbed in cold (Fig. 5 H–Q). The parameter estimates (SI Appendix, Table S6) suggest that in drought, the conversion rates of GA₅₃ to GA₄₄ and GA₁₉ to GA₂₀ are similar to those in control conditions, whereas the conversion rate of GA₄₄ to GA₁₉ is approximately doubled, providing an explanation as to why GA₁₉ concentrations are similar in drought and control conditions while GA₄₄ is lower in drought. The parameter estimates suggest that in cold the degradation rates of both GA₂₀ and GA₂₉ are higher than in control conditions, so that overall degradation is substantially increased in cold despite the GA2ox transcript levels being lower.

Thus, our modeling approach identifies specific GA enzyme activities impacted by drought and cold and explains the observed metabolite levels. Although in drought GA20ox transcript levels are lower than in control conditions (Fig. 5I), the reduced model predicts that GA20ox enzymes mediate the GA₄₄-to-GA₁₉ oxidation step at a higher rate. This prediction suggests that GA₁ synthesis in drought is similar to that in control conditions and that the lower GA₁ concentrations are caused by increased degradation mediated by increased GA2ox transcript levels. In cold, the GA2ox transcripts are expressed at lower levels than in control conditions, but the rate constant associated with GA2ox-mediated GA₂₀ degradation is increased, explaining why GA₁ concentrations are similar in control and cold conditions.

To test the surprising model prediction that enzyme activities (relative to their transcript levels) increase in response to cold and drought, we performed further experiments to measure the enzyme activity directly. We considered six reactions for each condition: GA₅₃, GA₄₄, and GA₁₉ were used to determine GA20ox activities, and GA₁, GA₂₉, and GA₈ were used to obtain GA2ox activities. The enzymes were extracted from 10 mm segments from the maize leaf growth zone. These data (SI Appendix, Fig. S19) were used to calculate the oxidation rate constants (SI Appendix contains details).

In agreement with the model predictions, the oxidation rate mediated by GA2ox increased in cold conditions (Fig. 6). For drought, there was a trend for higher oxidation rates associated with GA20ox activity, also in support of the model predictions. Therefore, these in vivo activity measurements support the in silico prediction of increased specific enzyme activities in cold (GA2ox) and drought (GA20ox) conditions.

Full details of the mathematical model are provided in the SI Appendix. We defined a cell-based model integrating growth, metabolism, and dilution (SI Appendix, section 2) and used this to derive a reduced model (SI Appendix, section 3). We simulated the reduced model by specifying growth using experimental measurements (SI Appendix, section 2.3.1) and used the metabolite and transcript data to estimate the reduced model parameters using Matlab’s lsqnonlin optimization algorithm (SI Appendix, section 2.3.2–4). All code and data are provided in a GitLab repository.

Plant material and growth conditions were as described in ref. (31) (cold) and ref. (30) (drought).

For hormone profiling, sampling, extraction, purification, and hormone metabolic profiling were performed as described in ref. (9). Data for the UBI::GA20OX-1 experiment (Fig. 4, Fig S5) were reprinted from ref (9 Current Biology, Vol: 22, Hilde Nelissen, Bart Rymen, Yusuke Jikumaru, Kirin Demuynck, Mieke Van Lijsebettens, Yuji Kamiya, Dirk Inze, and Gerrit T.S. Beemster, A Local Maximum in Gibberellin Levels Regulates Maize Leaf Growth by Spatial Control of Cell Division, p1183-1187, 2012, with permission from Elsevier.

Enzyme transcript levels were measured as described in ref. (9). Transcript data for the UBI::GA20OX-1 experiment (Fig. 4) are as published in ref. (9). In the drought–cold experiment (Fig. 5), GA20ox levels are the summation of the measured levels of GA20ox1, GA20ox2.1, and GA20ox2.2; GA3ox is the level of GA3ox2; and GA2ox levels are the summation of GA2ox3.1, GA2ox3.2, GA2ox4, GA2ox6.2, GA2ox7.1, and GA2ox7.3, each measurement being the mean of n = 3 to 7 replicates for positions 0 to 30, located around the GA maximum; 95 mm, located at the end of the growth zone; and of n = 1 to 3 replicates for positions 35 to 85 mm, where GA levels are relatively stable.

Kinematic analysis and measurements of division–zone length were performed as described in ref. (32). The cross-sectional areas were calculated by measuring the volumes of 1 cm segments of leaf.

Enzymes were extracted from 10 segments of the fourth maize leaf, and enzyme activities were measured directly as described in SI Appendix, section 1.