# Sensitivity analyses for improving sulfur management strategies in winter oilseed rape

Emilie Poisson, Sophie Brunel-Muguet, François Kauffmann, Jacques Trouverie, Jean-Christophe Avice, Alain Mollier, Sara Amancio

► To cite this version:

Emilie Poisson, Sophie Brunel-Muguet, François Kauffmann, Jacques Trouverie, Jean-Christophe Avice, et al.. Sensitivity analyses for improving sulfur management strategies in winter oilseed rape. PLoS ONE, 2018, 13 (9), pp.e0204376. 10.1371/journal.pone.0204376. hal-02184355

HAL Id: hal-02184355

https://normandie-univ.hal.science/hal-02184355

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Sensitivity analyses for improving sulfur management strategies in winter oilseed rape

Emilie Poisson<sup>1,*</sup>, Sophie Brunel-Muguet<sup>1,*</sup>*, François Kauffmann<sup>2,‡</sup>, Jacques Trouverie<sup>1,‡</sup>, Jean-Christophe Avice<sup>1,‡</sup>, Alain Mollier<sup>1,§</sup>

1 UMR 950 EVA (Ecophysiologie Végétale et Agronomie), Normandie Université, UNICAEN, INRA, SFR Normandie Végétal, Caen, France, 2 UMR 6139 Laboratoire de Mathématiques Nicolas Oresme, Normandie Université, UNICAEN, CNRS, Caen, France, 3 ISPA, Bordeaux Sciences Agro, INRA, Villenave d'Ormon, France

© These authors contributed equally to this work.‡ These authors also contributed equally to this work.
* sophie.brunel-muguet@unicaen.fr

<figure><img src="image_2.png" /></figure>

## Abstract

Because sulfur (S) depletion in soil results in seed yield losses and grain quality degradation, especially in high S-demanding crops such as oilseed rape (*Brassica napus* L.), monitoring S fertilisation has become a central issue. Crop models can be efficient tools to conduct virtual experiments under different fertilisation management strategies. Using the process-based model SuMoToRI, we aimed to analyse the impact of different S fertilisation strategies coupled with the variability observed in major plant characteristics in oilseed rape i.e. radiation use efficiency (RUE), carbon (C) allocation to the leaves (β) and specific leaf area (SLA) on plant performance-driven variables encompassing total biomass (TDW), S in the photosynthetic leaves (QS<sub>mobile.GL</sub>) and leaf area index (LAI<sub>GL</sub>). The contrasting S supply conditions differed in the amount of S (5 levels), and the timing of application (at bolting and/or at flowering, which included a fractioned condition). For this purpose, we performed a global sensitivity analysis (GSA) and calculated two sensitivity indices i.e. the Partial Raw Correlation Coefficient (PRCC) and the Sobol index. The results showed that whatever the timing of S supply, TDW, LAI<sub>GL</sub> and QS<sub>mobile.GL</sub> increased as S input increased. For a given S supply, there was no difference in TDW, LAI<sub>GL</sub> and QS<sub>mobile.GL</sub> between a single and a fractioned supply. Moreover, delaying the supply until flowering reduced the TDW and LAI<sub>GL</sub> whereas QS<sub>mobile.GL</sub> increased. Results showed that RUE had the greatest impact on TDW under all levels of S supply and all application timings, followed by β and SLA. RUE mostly impacted on QS<sub>mobile.GL</sub>, depending on S supply conditions, whereas it was the parameter with the least impact on LAI<sub>GL</sub>. Ultimately, our results provide strong evidence of optimised S fertilisation timings and plant characteristics that will guide producers in their agricultural practices by using specific varieties under constrained S fertilisation strategies.

<figure><img src="image_3.png" /></figure>

### OPEN ACCESS

**Citation:** Poisson E, Brunel-Muguet S, Kauffmann F, Trouverie J, Avice J-C, Mollier A (2018) Sensitivity analyses for improving sulfur management strategies in winter oilseed rape. PLoS ONE 13(9): e0204376. https://doi.org/10.1371/journal.pone.0204376

**Editor:** Sara Amancio, Universidade de Lisboa Instituto Superior de Agronomia, PORTUGAL

**Received:** June 19, 2018

**Accepted:** September 6, 2018

**Published:** September 20, 2018

**Copyright:** © 2018 Poisson et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files.

**Funding:** This work was performed, in partnership with the SAS PIVERT, within the framework of the French Institute for Energy Transition (Institut pour la Transition Energétique (ITE) Picardie Innovation Végétales, Enseignements et Recherches Technologiques (P.I.V.E.R.T.) (www.institut-pivert.com) selected as an Investment for the Future ("Investissements d'Avenir"). This work was

## Introduction

Crop models can be efficient tools to simulate and analyse agricultural practices, such as the impact of fertiliser management strategies on crop performances [1–6]. They provide a
supported, as part of the Investments for the Future by the French Government under the reference, ANR-001-01. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

framework to untangle interactions between plants and their environment, and often deliver quantitative outputs indicating thresholds to target specific objectives. In the current context of sustainable agriculture, which implies reduction of fertiliser inputs, modelling experiments are less expensive and easier to set up than field or greenhouse experiments but allow fertilisation constraints to be tested. In addition, most crop models have been used extensively to analyse the impacts of on-going global changes and to project crop yields worldwide [7–11].

Modelling approaches to assist crop management in low input systems can target fertilisation strategies (amount, timing and fractioning) or plant characteristics, which in both cases aim to alleviate the negative impact of decreased resources. The advantage of process-based models is that (i) they clarify plant characteristics and central processes such as photosynthesis, nutrient uptake, assimilation or remobilisation by means of equations and parameters, and (ii) they take into account the variability of plant features. Therefore, exploring the effects of the variation in plant parameters on model outputs can be achieved with sensitivity analyses (SAs) [12]. SAs are usually carried out to assess the sensitivity of model outputs with respect to the inputs i.e. parameters and/or variables (due to their estimation uncertainty or inherent variability). SAs with process-based models have several objectives: (i) to explored model functioning; (ii) to improve model calibration by adding or deleting processes (i.e. parameters and equations) which define the level of complexity in the model according to their impacts on output variations [13]; (iii) to quantitatively estimate the extent of the impacts of input variables (representing resource allocation for instance) on crop performance; and (iv) to identify the most influential plant parameters (representing plant features) according to their effect on output variations, and thus to help defining ideotypes (or virtual genotypes). Two categories of SAs are usually carried out. First, the local approach (LA) aims to study the variability of outputs induced by the small variation of a single parameter and/or variable (input) to the point when an estimate is chosen while keeping the other parameters fixed [14]. Second, a global sensitivity analysis (GSA) considers in contrast the whole variation range or domains of the inputs, which allows outputs to be analysed with regard to the combined variability [14]. GSA are commonly used with complex models where the conventional method (LA) would fail to represent the interactions between the numerous parameters and input variables [15].

In this study, we aimed to analyse the impact of different sulfur (S) fertilisation strategies coupled with the variability observed in major plant characteristics in oilseed rape. Like most of the Brassica species, oilseed rape is a high-S demanding crop with recommended inputs (RI) (Terres Inovia, http://www.terresinovia.fr/colza/cultiver-ducolza/fertilization/soufre/) of about 30 kg S.ha⁻¹ (corresponding to 75 kg SO₃.ha⁻¹) [16]. Sulfur is present in a wide range of metabolites such as methionine (essential amino acids (AA)) and cysteine (non-essential AA), proteins, glutathione and glucosinolates, which makes it an essential element for growth, development and resistance to abiotic and biotic stresses. The importance of tightly monitoring S fertilisation, especially in Brassica species, has emerged over recent decades when drastic environmental policies aimed at lowering industrial sulfur dioxide (SO₂) emissions [17] and in fine led to an increased risk of soil S oligotrophy. Several studies in oilseed rape have indicated that low S availability can drastically impact seed yield and corresponding components (e.g. thousand seed weight, pod number, seeds per pod), seed nutritional characteristics (e.g. oil and protein contents, and fatty acid (FA) and seed storage protein (SSP) compositions) as well as physiological quality of seeds (e.g. germination rates, storage capacity) [18–20]. In this context, modelling S requirements and allocation within the plant could be a helpful approach to correct deficiencies occurring during vegetative growth and the transition between vegetative and reproductive stages which are crucial phases for S storage and thus for providing S to growing reproductive sinks [21–23]. The process-based model SuMoToRI (Sulphur Model Towards Rapeseed Improvement) [2] simulates the dynamics of crop growth and S status until
the onset of pod formation. Its original feature is based on the prediction of an S-mineral pool that is assumed to remobilise to provide S requirements for growing sinks throughout the reproductive stage. The model is run with a limited number of plant parameters that refer to the potential leaf growth and carbon (C) demand by the leaves driven by climatic variables (temperature and radiation), genetic characteristics, C assimilation and S uptake as well as S allocation (among compartments) and partitioning (S-organic vs. S-mineral compounds).

The objective of this modelling analysis is to highlight the most suitable combinations of plant model parameters under contrasting fertilisation strategies which allowed high plant performances (biomass, leaf area or S content in leaves). For this purpose we performed a global sensitivity analysis (GSA) and determined two sensitivity indices i.e. the Partial Raw Correlation Coefficient (PRCC) and the Sobol index. Both indices are based on measures of importance (quantitative sensitivity indices) in contrast with coarse screening methods e.g. the Morris method [24], which is usually used to explore the behaviour of model outputs by varying a broad number of parameters, thus helping to detect the least influential ones. Because we aimed to focus on specific plant parameters, we used methods based on linear regression (i.e. PRCC) and on the decomposition of functional variance using Sobol indices. This allowed (i) estimation of the strength of linear associations between outputs and each input parameter after removing the linear effect of the other parameter (PRCC) and (ii) distinguishing and quantifying the principal and interaction effects of the parameter on the output variation using Sobol indices [25].

In our study, the GSA was performed under contrasting S supply conditions that differed in S availability (5 levels), date (at bolting or at flowering) and fractioning (once or twice during growth from the end of winter) of S inputs. Three plant parameters of the model SuMoToRI were targeted because they were shown to display variability in response to S [2] and/or different conditions that have been reported in the literature [26]. These parameters were the Radiation Use Efficiency (RUE, g DW.MJ⁻¹), the Specific Leaf Area (SLA, m².g DW⁻¹) and the C-leaf allocation coefficient (β, dimensionless). The outcomes for the GSAs are twofold (i) ranking the targeted plant parameters according to their impacts on plant performances, that are the total biomass production, the S in the leaves available for remobilisation towards growing pods (used as a proxy of seed quality) and the leaf area index (the central variable in the model) under different S supply strategies and (ii) identifying the most suitable combinations of plant parameter values (namely ideotypes) under these contrasting S supply strategies. This will eventually help producers to adjust agricultural practices by using specific varieties under constrained S fertilisation strategies.

# Materials and methods

## Overview of the SuMoToRI model

SuMoToRI was described extensively in Brunel-Muguet *et al.* [2]. Briefly, this process-based model predicts with daily time increments, the dynamics of the leaf area index (LAI<sub>GL</sub> which is the leaf area, LA<sub>GL</sub>, multiplied by plant density), the biomass, the S amounts and the fractions of organic and mineral S for three main plant compartments considered, namely the photo-synthetic leaves, which are simplified as a single Big Leaf (BL), the fallen leaves (FL) and the rest of the plant. The model considers the three environmental factors of temperature, Photo-synthetically Active Radiation (PAR, MJ m⁻²) and the amount of S taken up by the plant. The simulation period covers the end of vernalisation until the onset of pod formation. The mineral fraction in the leaves is estimated by the sulfate amount and this is used as an indicator of the potential for S remobilisation towards growing pods. The model is run with 23 plant
parameters, most with generic values, which describe potential leaf expansion, C-assimilation, allocation of C and S among the three compartments and S-partitioning (mineral vs. organic).

## Simulations and sensitivity analysis procedures

The model was used to predict plant growth and S status under several S supply conditions that were expected to highlight contrasting plant behaviour simulations. Then the GSA was performed (i) to rank three plant parameters according to their respective impact on plant performance and (ii) to determine the most suitable plant parameter value combinations under these special conditions. Our underlying questions were the following: (i) what are the most impacted outputs resulting from the variations in the targeted parameters? (ii) To what extent do S fertilisation conditions modulate the impact of the parameter variations on outputs? The GA procedure consists of the following steps:

**Step 1. Choice of targeted plant parameters and setting of their variation range.** The GSA was performed on three of the 23 parameters in order to assess their impact on biomass and S-content: the Radiation Use Efficiency (RUE), the Specific Leaf Area (SLA), and the β parameter, which indicates the allocation of C-assimilates to the leaves. These parameters were chosen (i) because value variations were observed in response to S availability for RUE (variability observed in the model calibration and evaluation datasets) and for the Specific Leaf Area (SLA), and the β coefficient (with variability for the validation dataset only) (unpublished data and [2]) and (ii) according to the literature and prior experiments including those reported from other species [26–28]. Following these observations, the mean, the standard deviation, the distribution profile and the truncation thresholds (minimum and maximum) were determined for each parameter (Table 1). All three parameter distributions were assumed to be uniform (not negative). All other model parameters were assumed constant (Table 1).

**Step 2. Calculation of S uptake functions.** The model was initially calibrated and evaluated with the cultivar Yudal, under two contrasting S supply conditions (High S, HS and Low S, LS) [29]. They corresponded respectively to 300 and 20 units (U) of sulfur trioxide (SO₃, kg·ha⁻¹), which were provided throughout the crop cycle following the relative addition rate nutrient-dosing system [30,31]. Consequently, measured plant S uptake dynamics for both extreme conditions were fitted to the Hill’s model (Eq 1):

$$QS = \frac{QS_{max} \times TT^n}{K_a + TT^n} + QS_{ini} \qquad (Eq 1)$$

where QSmax (mg S.plant⁻¹), K_a (°Cd) and n parameters were the three plant parameters describing the uptake process, QSini (mg S.plant⁻¹) the initial amount taken up by the plant at the end of vernalisation and TT, the thermal time (in °Cd) [32]. Three other intermediate S fertilisation levels were selected as follows (i) the recommended inputs (RI) for the whole crop cycle, matching 75U SO₃, (ii) 50U SO₃ (2/3 of the RI) and (iii) 37.5U SO₃ (half of the RI). These amounts were supplied once at the end of winter (GS30, bolting) and 20 days later (200 °Cd, GS60, flowering) [33] (Fig 1A), or split at two times in equal amounts at GS30 and GS60 (Fig 1B) according to the following Hill’s model (Eq 2):

$$QS = \frac{QS_{max1} \times TT^{n1}}{K_{a1} + TT^{n1}} + \frac{QS_{max2} \times TT^{n2}}{K_{a2} + TT^{n2}} + QS_{ini} \qquad (Eq 2)$$

The Hill’s model parameters were determined for the 15 S-supply conditions i.e. 5 amounts x 3 timings (once at GS30, once at GS60 and twice at GS30 and GS60 for the fractioned condition) (Table 2). Initial values were the same for the 15 conditions, assuming that the plants
<table><caption>Table 1. Model parameters and initial values.</caption><thead><tr><th>Symbol</th><th>Definition</th><th>Value</th><th>Unit</th><th>Source</th></tr></thead><tbody><tr><td>sd</td><td>Sowing density</td><td>40</td><td>Plant.m<sup>-2</sup></td><td></td></tr><tr><td colspan="5"><strong>PAR interception</strong></td></tr><tr><td>k</td><td>PAR extinction coefficient</td><td>0.75</td><td>m<sup>2</sup> m<sup>-2</sup></td><td>Bonhomme et al. 1982</td></tr><tr><td colspan="5"><strong>Potential leaf growth</strong></td></tr><tr><td>LA<sub>0</sub></td><td>Initial leaf area of photosynthetic leaves</td><td>0.015</td><td>m<sup>2</sup> plant<sup>-1</sup></td><td>Estimated</td></tr><tr><td>LA<sub>max</sub></td><td>Leaf area expansion parameters</td><td>0.2</td><td>m<sup>2</sup> plant<sup>-1</sup></td><td></td></tr><tr><td>K</td><td></td><td>872.96</td><td>°Cd<sup>-1</sup></td><td></td></tr><tr><td>N</td><td></td><td>6.31</td><td>dimensionless</td><td></td></tr><tr><td colspan="5"><strong>C acquisition and plant offer</strong></td></tr><tr><td>PARabs<sub>ini</sub></td><td>Initial absorbed PAR</td><td>0</td><td>MJ m<sup>-2</sup></td><td>Estimated</td></tr><tr><td>TDW<sub>ini</sub></td><td>Initial total dry weight</td><td>0.576</td><td>g DW plant<sup>-1</sup></td><td></td></tr><tr><td>DW<sub>FL ini</sub></td><td>Initial dry weight of fallen leaves</td><td>0</td><td>g DW plant<sup>-1</sup></td><td></td></tr><tr><td><strong>RUE</strong></td><td><strong>Radiation use efficiency</strong></td><td><strong>1.6–4.6</strong></td><td><strong>g DW MJ<sup>-1</sup></strong></td><td></td></tr><tr><td>aLDW<sub>FL</sub></td><td>Parameters of the function describing the time progression of LDW<sub>FL</sub></td><td>0.0092</td><td>g DW plant<sup>-1</sup> °Cd<sup>-1</sup></td><td></td></tr><tr><td>bLDW<sub>FL</sub></td><td></td><td>0.0043</td><td>dimensionless</td><td></td></tr><tr><td colspan="5"><strong>C allocation to leaves</strong></td></tr><tr><td>β</td><td>Coefficient of DW allocation to the leaves</td><td>0.10–0.72</td><td>dimensionless</td><td>Estimated</td></tr><tr><td colspan="5"><strong>C demand of green leaves</strong></td></tr><tr><td>LDW<sub>GL ini</sub></td><td>Initial dry weight of green leaves</td><td>0.448</td><td>g DW plant<sup>-1</sup></td><td>Estimated</td></tr><tr><td><strong>SLA</strong></td><td><strong>Specific leaf area</strong></td><td><strong>0.008–0.034</strong></td><td><strong>m<sup>2</sup> g DW<sup>-1</sup></strong></td><td></td></tr><tr><td colspan="5"><strong>Growth S demand</strong></td></tr><tr><td>αGL</td><td>Parameters to estimate critical S content in GL as a function of LDWGL</td><td>5.11</td><td>mg S plant<sup>-1</sup></td><td>Estimated</td></tr><tr><td>βGL</td><td></td><td>-0.52</td><td>dimensionless</td><td></td></tr><tr><td>αrest</td><td>Parameters to estimate critical S content in the rest of the plant as a function of DWrest</td><td>1.83</td><td>mg S plant<sup>-1</sup></td><td></td></tr><tr><td>βrest</td><td></td><td>-0.004</td><td>dimensionless</td><td></td></tr><tr><td colspan="5"><strong>Potential mobile S allocation</strong></td></tr><tr><td>εpot</td><td>Coefficient of potential repartition of mobile S to the leaves</td><td>0.8</td><td>dimensionless</td><td>Estimated</td></tr><tr><td colspan="5"><strong>S uptake</strong></td></tr><tr><td>QS<sub>TOT ini</sub></td><td>Initial total S uptake</td><td>6.78</td><td>mg S plant<sup>-1</sup></td><td>Estimated</td></tr><tr><td>QS<sub>GL ini</sub></td><td>Initial S in green leaves</td><td>5.76</td><td>mg S plant<sup>-1</sup></td><td></td></tr><tr><td>QSrest<sub>ini</sub></td><td>Initial S in fallen leaves</td><td>1.017</td><td>mg S plant<sup>-1</sup></td><td></td></tr></tbody></table>

Parameters in bold are used for the GSAs (RUE, β and SLA). Values range between minimum and maximum.

https://doi.org/10.1371/journal.pone.0204376.t001

were previously grown under similar environmental conditions and optimal S nutritional supply (75 U SO<sub>3</sub>). Therefore, the initial amount of total S (QS<sub>ini</sub>) was determined relative to both HS and LS conditions by averaging the values obtained with the respective Hill's adjustments for HS and LS. Then linear regressions were used to calculate the initial total dry weight (TDW<sub>ini</sub>), the initial leaf dry weight of green leaves (LDW<sub>GLini</sub>) and the initial amount of S in green leaves (QS<sub>GLini</sub>) (Table 1).

**Step 3. Simulations of plant performances under contrasting S-supply conditions.**

Simulations were performed under 15 different S-supply conditions that differed in terms of S-amount and timing of application (including one fractioned condition) (Table 2). A single
<figure><img src="image_4.png" /><figcaption>**Fig 1.** Hills' kinetics of S uptake according to the 15 S-supply conditions. Fitted Hill's model of S amounts as a function of thermal time (°Cd) for (A) the single S supplies at GS30, (B) at GS60 and (C) the fractioned S supplies at GS30 and GS60. Tb: base temperature.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0204376.g001

set of parameter values was used except for the three targeted parameters (RUE, SLA and β), which were randomly selected within their respective distribution (Tables 1 and 2). A single climatic dataset for simulations and the GSAs was used, based on daily average temperature and radiation from 2005 to 2015 in Saint-Martin-de-Hinx (43°34'57"N, 1°16'10"W) (data from the CLIMATIK platform, https://intranet.inra.fr/climatik_v2, S1 Fig, supporting information). The dates of initialisation for the simulations were determined with two criteria i.e. temperature (10 consecutive days above 10°C) and day length (threshold of 11.35 hours), so as to closely match the conditions at the end of vernalisation to start running the model [29].

<table><caption>**Table 2.** Hill's model parameters for the 15 S-supply conditions.</caption><thead><tr><th>Total S rate supply (Units: kg.ha<sup>-1</sup>)</th><th>Definition</th><th>Q<sub>Smax1</sub></th><th>n<sub>1</sub></th><th>K<sub>a1</sub></th><th>Q<sub>Smax2</sub></th><th>n<sub>2</sub></th><th>K<sub>a2</sub></th></tr></thead><tbody><tr><td colspan="8"><strong>Supply at the end of winter</strong></td></tr><tr><td>20</td><td rowspan="5">Corresponding to a supply of S at the end of winter (GS30: beginning of simulation).</td><td>14.5</td><td>9.1</td><td>654</td><td>/</td><td>/</td><td>/</td></tr><tr><td>37.5</td><td>26.6</td><td>9.0</td><td>656</td><td>/</td><td>/</td><td>/</td></tr><tr><td>50</td><td>35.2</td><td>8.9</td><td>658</td><td>/</td><td>/</td><td>/</td></tr><tr><td>75</td><td>52.4</td><td>8.6</td><td>661</td><td>/</td><td>/</td><td>/</td></tr><tr><td>300</td><td>207.5</td><td>6.7</td><td>686</td><td>/</td><td>/</td><td>/</td></tr><tr><td colspan="8"><strong>Supply 20 days delayed</strong></td></tr><tr><td>20</td><td rowspan="5">Corresponding to a supply of S 20 days after the end of winter (GS60).</td><td>14.5</td><td>24.6</td><td>854</td><td>/</td><td>/</td><td>/</td></tr><tr><td>37.5</td><td>26.6</td><td>24.4</td><td>856</td><td>/</td><td>/</td><td>/</td></tr><tr><td>50</td><td>35.2</td><td>24.2</td><td>858</td><td>/</td><td>/</td><td>/</td></tr><tr><td>75</td><td>52.4</td><td>24.0</td><td>861</td><td>/</td><td>/</td><td>/</td></tr><tr><td>300</td><td>207.5</td><td>20.8</td><td>886</td><td>/</td><td>/</td><td>/</td></tr><tr><td colspan="8"><strong>Fractioned S application</strong></td></tr><tr><td>20</td><td rowspan="5">Corresponding to a fractioned supply of S at the end of winter and 20 days after the end of winter (GS30 +GS60).</td><td>7.3</td><td>18.3</td><td>327</td><td>7.3</td><td>18.3</td><td>654</td></tr><tr><td>37.5</td><td>13.3</td><td>18.0</td><td>328</td><td>13.3</td><td>18.0</td><td>656</td></tr><tr><td>50</td><td>17.6</td><td>17.7</td><td>329</td><td>17.6</td><td>17.7</td><td>658</td></tr><tr><td>75</td><td>26.2</td><td>17.3</td><td>330</td><td>26.2</td><td>17.3</td><td>660</td></tr><tr><td>300</td><td>103.7</td><td>13.3</td><td>342</td><td>103.7</td><td>13.3</td><td>685</td></tr></tbody></table>

Equations for S-uptake adjustments are given in the material and methods section.

https://doi.org/10.1371/journal.pone.0204376.t002
**Step 4. Computation of the sensitivity indices: partial rank correlation coefficient.** The partial rank correlation coefficient (PRCC) measures the strength of the linear associations between the output and each input parameter, after removing the linear effect of the other parameters. These rank-based measures are part of the so-called sampling-based global sensitivity analysis method. The PRCC varies between -1 and +1 and it quantifies the links between input factors and output variables as well as the direction of the relationship. We used a “Latin Hypercube Sample” (LHS) for generating the sample of parameter combinations [34], which allows the precision of the sensitivity indices to be increased. Overall 3,000 simulations (200 repetitions x 15 conditions) for each output were performed, thus allowing 200 combinations of parameters to be tested under the 15 S-supply conditions.

**Step 5. Computation of the sensitivity indices: functional decomposition of variance with Sobol indices.** The Sobol method is a variance-based method that uses a variance ratio to estimate the importance of parameters [35,36]. Two main sensitivity indices are defined for each parameter. The first order or main order of sensitivity index measures the average effect of one parameter on one model output, without taking into account the interaction effects with the other parameters. The second order expresses the sensitivity of the model to the interactions between the parameters. The sum of both indices provides the total effect index. Prior to their estimation, we used the Monte Carlo-based sampling method, which implied n (d+2) model evaluations (N) in which n = 2000 is the size of an initial Monte Carlo sample (number of repetitions) and d = 3 is the number of targeted parameters [13]. Therefore, the index estimations required 10⁴ model evaluations (N) per S supply condition and for each output. Overall, for each output variable, 10⁴ x 15 simulations were performed under the 15 S-supply conditions.

## Software

The model was run with R (version 3.4.1) [37] with additional R packages including *pse* (Latin Hypercubes), *sensitivity* (sensitivity analysis), *ggplot2*, *ggthemes* and *reshape2* (visualisation).

## Results

Three representative outputs were selected to illustrate the variations in the model inputs and these comprised the biomass (TDW), the leaf area index (LAI<sub>GL</sub>) and the S-remobilisation process (QS<sub>mobile,GL</sub>).

### Plant performances and range of variations at the onset of pod formation

Fig 2 represents the extent of variation in the 3 outputs (i.e. TDW, LAI<sub>GL</sub> and QS<sub>mobile,GL</sub>) at the onset of pod formation for the 15 S supply conditions according to random draws of 200 combinations of the 3 parameters tested (i.e. RUE, SLA and β).

Whatever the timing of the S supply (GS30, GS60 and GS30+60, which is the fractioned condition), the three outputs (TDW, LAI<sub>GL</sub> and QS<sub>mobile,GL</sub>) increased with a higher S supply but for TDW the increase was lower (Fig 2). The median TDW ranged from 32.9–39.7 g.plant<sup>-1</sup> with extreme values observed for 20U and 300U, respectively, when provided at GS30 (bolting). The S supply rate effect was stronger on LAI<sub>GL</sub> when provided once at GS30 and twice at GS30 +60 (with values being more than twice as high with 300U compared to 20U) and on QS<sub>mobile,GL</sub>, irrespective of the timing of the S supply (with the expected increase being up to 10<sup>4</sup> times higher with 300U than with 20U).

For a given S supply rate there was a slight increase in the TDW, LAI<sub>GL</sub> and QS<sub>mobile,GL</sub> (Fig 2A, 2B and 2C) when supplying S twice at GS30+60 compared to once at GS30. However, delaying S fertilisation until GS60 noticeably reduced the TDW and LAI<sub>GL</sub> (Fig 2A and 2B),
<figure><img src="image_5.png" /><figcaption>**Fig 2.** Plant performances and range of variations at the onset of pod formation for the 15 S-supply conditions according to random draws of 200 combinations of the 3 parameters tested. Variation of the TDW (g.plant⁻¹) (A), LAI<sub>GL</sub> (m².m⁻²) (B) and QS<sub>mobile,GL</sub> (mg S.plant⁻¹) (C) for the 15 S-amounts x timing conditions, obtained from simulations with 200 combinations of the 3 parameters tested (RUE, β and SLA). Data are represented by a box plot and through their quartiles; the bottom and top of the box are the first or lower (Q1) and third or upper quartiles (Q3), the band inside the box is the second quartile (the median, Q2) and maximum and minimum values are at the end of the vertical bars. Extreme data are represented by point.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0204376.g002

regardless of the S amount, and contrasted with the increase in QS<sub>mobile,GL</sub> (Fig 2C). The effect of the S amount when provided at GS60 (flowering) had less of an impact on TDW (25.8 to 27.3 g.plant⁻¹) and LAI<sub>GL</sub> (1.5 to 2.5 m².m⁻²) than supplying S at GS30 (32.9 to 39.7 for TDW and 2.3 to 6.6 for LAI<sub>GL</sub>) or at GS30+60 (35.6 to 41.4 for TDW and 2.3 to 6.8 for LAI<sub>GL</sub>).

The simultaneous variation in the 3 parameters generated a high range of variation in the three outputs (Fig 2A, 2B and 2C). For TDW and LAI<sub>GL</sub> the extent of variation resulting from combinations of the three parameters was more pronounced with the GS30 (bolting) and GS30+60 timings than in the GS60 timing (flowering) under all S supply rates, except in the case of LAI<sub>GL</sub> where lower variation was observed with lower S amounts. For instance, the extent of variation between the minimum (min) and the maximum (max) TDW values were *ca.* 41.7 and 27.8 respectively for 20U at GS30 or GS60 (Fig 2A). For LAI<sub>GL</sub>, the extent of variation between max and min values was about 4.2 and 5.9 respectively for 20U and 300U at GS30 and about 2.5 and 3.6 respectively for 20U and 300U at GS60 (Fig 2B). In contrast, the extent of variation for QS<sub>mobile,GL</sub> was much lower than for the other two parameters, but it
was intensified under GS60 timing conditions with a similar trend for each of the five S amounts, whereas under the GS30 and GS30+60 timing conditions, the extent of variation for QS<sub>mobile,GL</sub> only increased with 75U and 300U (Fig 2C).

**Impact of variation in plant parameters on plant performance and S status**

To investigate the specific impact of the 3 parameters (RUE, SLA and β), global sensitivity analyses were performed. They allowed sensitivity indices (PRCC and Sobol indices) to be calculated for the three targeted outputs under the 15 conditions. Furthermore, the Sobol indices indicated and quantified the interactions between the three parameters. By doing so, we sought to rank the plant parameters according to the impact of their variations on plant performance and to see whether the S amounts x timing designs could interfere in this ranking. The results are supported by outputs that illustrate the main processes included in the model.

**The variation in RUE was the main driver of TDW and displayed high stability across the conditions.** Plant biomass at the onset of pod formation was strongly positively impacted by the variation in RUE when PRCC values were close to 1 according to the S amounts and timing conditions (Fig 3A). Furthermore, the increase in β and SLA was positively correlated with the increase in TDW when PRCC values ranged from 0.71 to 0.86 and from 0.72 to 0.86 for β and SLA respectively under all conditions. The small range of variations in the PRCC values of RUE, β and SLA indicated that S-amounts and fertilisation timing conditions did not influence the respective impact of these three parameters on TDW. Similar conclusions could be drawn with the Sobol indices calculated for the three parameters (Fig 3B). Among the parameters, RUE had the most impact with a mean main index of 0.85 across all conditions together, which was in contrast to the much lower main index values for β and SLA. The interaction index values were low for the three parameters indicating no interaction was detected between them (Fig 3B).

**LAI<sub>GL</sub> was mainly impacted by the C-allocation to the leaves.** Both sensitivity indices indicated that mainly SLA and β highly influenced LAI<sub>GL</sub> with different patterns depending on the levels of S fertilisation and timing conditions (Fig 4A and 4B). Regarding the PRCC indices, the mean index value of RUE reached 0.52 across all conditions together and decreased in the fractioned condition (GS30+60) with S inputs below 300U. The positive β index was the highest and was unaffected by the S amount or timing conditions. Finally, SLA index values were variable under changing S amounts in the GS30 and GS30+60 conditions, with a decrease in the index under higher S inputs. The Sobol indices confirmed that the influence of RUE variation on LAI<sub>GL</sub> was the lowest, irrespective of the S amount and the timing conditions (Fig 4B). In contrast, the variation in β and SLA led to contrasting responses to the S amount and timing conditions. Under GS30 and GS30+60 conditions, the higher the S amount, the higher the impact of the variation in β on LAI<sub>GL</sub>. For instance, the Sobol total index under the GS30 conditions increased from 0.3 to 0.7 as the amount of S increased. In contrast, under GS30 and GS30+60 conditions, the higher the S amount, the lower the impact of SLA on LAI<sub>GL</sub> whose Sobol total index decreased from 0.6 to 0.3 under the GS30 conditions. Under GS60 conditions, no effect of increased S amounts was observed on the sensitivity of LAI<sub>GL</sub> to the variation in β and SLA. Interaction indices were low for the three parameters meaning that there were no tight interactions between them, regardless of the S amounts and timing conditions (Fig 4B).

**QS<sub>mobile,GL</sub> was strongly impacted by RUE and the impact of β and SLA were depending on S supply.** The impact of the three parameters on QS<sub>mobile,GL</sub> was very variable according to the S amounts and the timing conditions (Fig 5A and 5B). Both the PRCC and Sobol indices indicated the strong influence of the variation in RUE on QS<sub>mobile,GL</sub>. The PRCC values highlighted similar patterns of sensitivity of QS<sub>mobile,GL</sub> to variations in RUE, β and SLA
<figure><img src="image_6.png" /><figcaption>**Fig 3.** Impact of variation in plant parameters on TDW described by PRCC and Sobol indices. Impact of RUE, β and SLA calculated by PRCC (A) and by Sobol indices (B) on TDW for the 15 S-amounts x timing conditions. For Sobol indices, total indice ±SE, main indice ±SE and interaction (= total-main) are presented for each parameter and each condition (n = 15).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0204376.g003

depending on the S amounts and timing conditions (Fig 5A). PRCC values were negative and tended to get closer to 0 with low S amounts under the GS30 and the GS30+60 conditions or with high amounts of S under GS60 conditions. Therefore, the higher the S amounts when
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Impact of variation in plant parameters on LAI<sub>GL</sub> described by PRCC and Sobol indices. Impact of RUE, β and SLA calculated by PRCC (A) and by Sobol indices (B) on LAI<sub>GL</sub> for the 15 S-amounts x timing conditions. For Sobol indices, total indice ±SE, main indice ±SE and interaction (= total-main) are presented for each parameter and each condition (n = 15).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0204376.g004

supplied at GS30 or GS30+60, the higher the negative impact of the variation in RUE, β and SLA on QS<sub>mobile,GL</sub>. In contrast, under GS60 conditions the impacts of the variations in RUE, β and SLA was not so dependent on the S amounts because the PRCC values were more stable,
<figure><img src="image_8.png" /><figcaption>**Fig 5.** Impact of variation in plant parameters on QS<sub>mobile</sub> described by PRCC and Sobol indices. Impact of RUE, β and SLA calculated by PRCC (A) and by Sobol indices (B) on QS<sub>mobile,GL</sub> for the 15 S-amounts x timing conditions. For Sobol indices, total indice ±SE, main indice ±SE and interaction (= total-main) are presented for each parameter and each condition (n = 15).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0204376.g005

with mean values of *ca*. -0.79, -0.54 and -0.31 for RUE, β and SLA, respectively, with all condi-
tions combined. Regarding the Sobol indices, RUE was the most influential parameter in
terms of total and main effects (Fig 5B) for most of the S-amounts and timing conditions,
except for the following fertilisation designs: 300U at GS60 (in this case β was the most influential), and 37.5U and 50U at GS30+60 (in these cases RUE and β were equally influential). The total index for β decreased concomitantly with the increase in the S amount at GS30 and GS30+60 (Fig 5A and 5B) as also observed for SLA (Fig 5B), which contradicts the PRCC indices. However, the impact of the variation in the three parameters was less influenced by the S amount when provided at GS60 with its lower Sobol indices, except in the case of 300U. The interaction indices varied for the three parameters across all S amounts and timing conditions. The interaction values decreased as the S amounts increased under GS30 and GS30+60 conditions when the S amounts were higher than 20U, but the same trend was not seen at GS60 where lower and stable interaction indices across the S amounts were observed.

## Identification of optimised parameter combinations under different S-fertilisation strategies

The PRCC analyses were conducted from a random draw of 200 combinations of the 3 parameters associated with an output value. Therefore, we could identify the most relevant plant parameter combinations (i.e. optimised parameter combinations) that enabled the best plant performances under the different fertilisation strategies tested. The best plant performances were associated with high biomass (TDW) and green leaf area (LAI<sub>GL</sub>) to maximise plant development and light interception and high levels of mobile S (QS<sub>mobile,GL</sub>) to allow adequate S storage in pods. The condition where 75U of S was applied at GS30 was used as the S supply reference condition (RC). The best plant performances obtained under RC were associated to optimised parameter combinations (RC<sub>opt</sub>) which were thus compared to other optimised parameter combinations under the other S supply conditions (Table 3).

<table><caption>Table 3. Highest values of the three outputs (TDW, LAI<sub>GL</sub> and QS<sub>mobile,GL</sub>) for each of 15 S-supply conditions obtained with optimized parameters combinations.</caption><thead><tr><th rowspan="2"></th><th colspan="4">TDW</th><th colspan="4">LAI<sub>GL</sub></th><th colspan="4">QS<sub>mobile,GL</sub></th></tr><tr><th>Δ TDW</th><th>RUE</th><th>SLA</th><th>β</th><th>Δ LAI<sub>GL</sub></th><th>RUE</th><th>SLA</th><th>β</th><th>Δ QS<sub>mobile,GL</sub></th><th>RUE</th><th>SLA</th><th>β</th></tr></thead><tbody><tr><td>Total S rate supply (kg.ha<sup>-1</sup>)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>RC<sub>opt</sub> (75 U of SO<sub>3</sub> at GS30)</td><td>56.4</td><td>4.58</td><td>0.016</td><td>0.54</td><td>6.49</td><td>2.59</td><td>0.018</td><td>0.65</td><td>25.68</td><td>1.61</td><td>0.022</td><td>0.11</td></tr><tr><td>Supply at GS30 (bolting)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>20</td><td>- 4%</td><td>0%</td><td>+ 83%</td><td>- 19%</td><td>- 24%</td><td>- 35%</td><td>+ 79%</td><td>- 26%</td><td>- 89%</td><td>=</td><td>- 61%</td><td>+ 17%</td></tr><tr><td>37.5</td><td>- 2%</td><td>- 3%</td><td>+ 91%</td><td>- 15%</td><td>- 6%</td><td>- 33%</td><td>+ 79%</td><td>- 3%</td><td>- 71%</td><td>+ 13%</td><td>- 53%</td><td>+ 49%</td></tr><tr><td>50</td><td>=</td><td>- 1%</td><td>+ 53%</td><td>+ 26%</td><td>- 2%</td><td>- 20%</td><td>+ 8%</td><td>- 1%</td><td>- 48%</td><td>=</td><td>- 45%</td><td>+ 83%</td></tr><tr><td>300</td><td>+ 13%</td><td>- 1%</td><td>+ 83%</td><td>+ 5%</td><td>+ 8%</td><td>+ 55%</td><td>- 43%</td><td>+ 3%</td><td>+ 413%</td><td>+ 2%</td><td>- 49%</td><td>+ 60%</td></tr><tr><td>Supply at GS60 (flowering)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>20</td><td>- 28%</td><td>- 4%</td><td>+ 110%</td><td>+ 25%</td><td>- 53%</td><td>- 31%</td><td>+ 87%</td><td>- 14%</td><td>- 75%</td><td>+ 8%</td><td>- 19%</td><td>+ 127%</td></tr><tr><td>37.5</td><td>- 27%</td><td>0%</td><td>+ 65%</td><td>- 33%</td><td>- 50%</td><td>- 21%</td><td>+ 27%</td><td>+ 8%</td><td>- 42%</td><td>+ 5%</td><td>- 51%</td><td>+ 107%</td></tr><tr><td>50</td><td>- 28%</td><td>- 2%</td><td>+ 2%</td><td>- 10%</td><td>- 46%</td><td>- 31%</td><td>+ 79%</td><td>+ 4%</td><td>- 19%</td><td>+ 7%</td><td>+ 19%</td><td>+ 95%</td></tr><tr><td>75</td><td>- 24%</td><td>0%</td><td>+ 97%</td><td>+ 11%</td><td>- 44%</td><td>- 33%</td><td>+ 71%</td><td>- 16%</td><td>+ 31%</td><td>+ 3%</td><td>- 58%</td><td>+ 181%</td></tr><tr><td>300</td><td>- 20%</td><td>- 1%</td><td>+ 97%</td><td>+ 3%</td><td>- 35%</td><td>- 27%</td><td>+ 81%</td><td>- 21%</td><td>+ 436%</td><td>+ 4%</td><td>- 61%</td><td>+ 130%</td></tr><tr><td>Fractioned S supply GS30+60</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>20</td><td>+ 6%</td><td>- 2%</td><td>+ 97%</td><td>+ 4%</td><td>- 37%</td><td>- 1%</td><td>+ 90%</td><td>+ 4%</td><td>- 98%</td><td>+ 97%</td><td>- 63%</td><td>+ 178%</td></tr><tr><td>37.5</td><td>+ 10%</td><td>0%</td><td>+ 41%</td><td>+ 20%</td><td>- 7%</td><td>- 35%</td><td>+ 42%</td><td>+ 5%</td><td>- 77%</td><td>+ 5%</td><td>- 63%</td><td>+ 46%</td></tr><tr><td>50</td><td>+ 12%</td><td>- 2%</td><td>+ 69%</td><td>- 44%</td><td>+ 2%</td><td>- 33%</td><td>- 2%</td><td>- 5%</td><td>- 41%</td><td>+ 6%</td><td>- 60%</td><td>+ 9%</td></tr><tr><td>75</td><td>+ 12%</td><td>- 1%</td><td>- 17%</td><td>+ 7%</td><td>+ 9%</td><td>- 13%</td><td>+ 2%</td><td>+ 1%</td><td>+ 7%</td><td>+ 2%</td><td>- 43%</td><td>+ 17%</td></tr><tr><td>300</td><td>+ 18%</td><td>0%</td><td>+ 10%</td><td>+ 24%</td><td>+ 13%</td><td>+ 51%</td><td>- 53%</td><td>+ 4%</td><td>+ 468%</td><td>+ 6%</td><td>- 60%</td><td>+ 32%</td></tr></tbody></table>

Figures in bold give the values of the outputs and parameters under the reference condition (RC) with optimized combinations (RC<sub>opt</sub>) which corresponds to 75 U of SO<sub>3</sub> at GS30. The other figures are given as the proportion of increase or decrease compared to the values under RC<sub>opt</sub>. The equal symbol (=) means there was no difference in values with RC<sub>opt</sub>.

https://doi.org/10.1371/journal.pone.0204376.t003
**Early and fractioned S supplies allowed the highest biomass production but different optimised plant parameter combinations.** When applied once at GS30, lower amounts of S than the reference condition (<75 U) led to small TDW variations even when the parameter combinations were optimised (-4% to -2% of RC<sub>opt</sub>). A higher S supply rate than in the RC<sub>opt</sub> (>75 U) resulted in a 13% increase in TDW with optimisation of the parameter combinations (Table 3). Whatever the S supply rate (sub and supra), the RUE values remained almost unchanged to those under RC<sub>opt</sub> with a maximum threshold value of 4.58 MJ.m<sup>-2</sup> under RC<sub>opt</sub>. In contrast, the SLA and β values were drivers of TDW variation because their variations were significant, ranging from +53% to +91% for SLA, irrespective of the amount of S and from -19% to +26% for β with increases only observed when S amount were above the RC<sub>opt</sub>. Under 50 U, the TDW remained unchanged if the SLA increased by 53% to RC<sub>opt</sub> and β increased by 25% compared to RC<sub>opt</sub> (75 U at GS30). When applied at GS60, there was a decrease in TDW from -20% to -28% compared to RC<sub>opt</sub> under all levels of S supply. When S was fractioned (GS30+60), the TDW increased from 6% to 18% under the extreme S supplies (20 and 300 U respectively) compared to RC<sub>opt</sub> (75U at GS30). The fractioned condition led to the highest TDW increase even under sub amounts (50U) if SLA increased by 69% and β values were almost two times lower than under RC<sub>opt</sub> (Table 3).

**The fractioned S supply allowed substantial increases in LAI<sub>GL</sub> with moderate changes in the plant parameters.** When S was supplied at GS30, LAI<sub>GL</sub> decreased by 24% to 2% compared to RC<sub>opt</sub> for S amounts below RC (<75 U) and increased by 8% with 300 U. The SLA values were higher whereas the RUE and β values were lower than under the RC<sub>opt</sub>, except under 300 U where SLA decreased by 43%, unlike RUE and β, which increased by 55% and 3%, respectively. Whatever the S amount, β was the least impacting driver with the smallest variations (-26% to 3%) while the ranges of values for RUE and SLA were much wider, with -35% to 55% for RUE and -43% to 79% for SLA (Table 3). For S supplies at GS60, despite optimised parameter combinations, LAI<sub>GL</sub> decreased regardless of the S amount, but as shown in Fig 2B the higher the S amount, the lower the decrease. The RUE values also decreased at GS60 and the SLA values increased compared to their values under RC<sub>opt</sub> (from -33% to -21% for RUE and from 27% to 87% for SLA, Table 3).

Under the fractioned supply (GS30+60), LAI<sub>GL</sub> increased with 50U, 75U and 300U by 2%, 9% and 13%, respectively, with optimised parameters combinations. For S amounts lower than RC (<75U), the optimised parameter combinations led to a reduction in LAI<sub>GL</sub> (or an insignificant increase, as observed under 50U).

**High QS<sub>mobile,GL</sub> was driven by the S amount, whatever the timing of supply, with little impact on optimised parameter combinations.** Whatever the timing of supply (at GS30, GS60 or GS30+60), QS<sub>mobile,GL</sub> decreased for S amounts lower than the RC (<75 U) and increased when the S supply was higher than RC. As expected, the increase was significant for the 300 U supplies at GS30+60 (up to 468% higher than under RC<sub>opt</sub>). The RUE and β values increased regardless of the S supply and timing compared to RC, while the SLA values decreased, except at GS60 with 50U. However, RUE values displayed the lowest increases compared to SLA and β for a given S condition. Interestingly, for amounts lower than RC (<75U), the increase in QS<sub>mobile,GL</sub> was higher when S was supplied once at GS60 (+31%) than twice under GS30+60 (+7%).

# Discussion

## Plant parameter rankings according to their impacts on plant performance under different S-fertilisation strategies

The three parameters (RUE, SLA and β) selected for the SAs are associated with C-metabolism and functioning. They were shown to be impacted by S supply [2], thus leading to conclude to
<table><caption>Table 4. Ranking of the 3 parameters according to their respective impact on TDW, LAI<sub>GL</sub> and QS<sub>mobile,GL</sub> according to S-supply condition.</caption><thead><tr><th>Output variable</th><th>Parameter ranking according to its impact</th><th>S-fertilization strategy</th><th>Correlation</th><th>Dynamics</th></tr></thead><tbody><tr><td>TDW</td><td>RUE > β ≈ SLA</td><td>whatever S amount and timing of application</td><td>RUE, SLA and β positively correlated</td><td>• Stable impact of increased RUE whatever S amount and timing<br>• More stable impact of increased SLA and β for supplies at GS60</td></tr><tr><td rowspan="2">LAI<sub>GL</sub></td><td>SLA > β > RUE</td><td>with 20U at GS30 and at GS30+60</td><td rowspan="2">RUE, SLA and β positively correlated</td><td rowspan="2">• Stable impact of increased β whatever S amount and timing<br>• Lower impact of increased SLA with higher S rate at GS30 and GS30+60<br>• Slightly higher impact of increased RUE with S rate at GS30 and lower impact of increased RUE with S rate up to 75U at GS30+60</td></tr><tr><td>β > SLA > RUE</td><td>other conditions</td></tr><tr><td rowspan="3">QS<sub>mobile,GL</sub></td><td>β > RUE > SLA</td><td>with 37.5U at GS30 with 300U at GS60</td><td rowspan="3">RUE, SLA and β negatively correlated<br>except RUE and β with 20U at GS30 and with 50U at GS30+60</td><td rowspan="3">• Higher impact of increased RUE, SLA and β with S rates especially at GS30 and GS30+60</td></tr><tr><td>RUE > SLA > β</td><td>with 20U and 300U at GS30</td></tr><tr><td>RUE > β > SLA</td><td>other conditions</td></tr></tbody></table>

This synthesis including results of the Sobol (main indices) and PRCC analyses.

https://doi.org/10.1371/journal.pone.0204376.t004

interactions between the C and S-related processes within the model. Our analysis highlighted that impacts on the variations (range and direction) of the outputs (TDW, LAI<sub>GL</sub> and QS<sub>mobile,GL</sub>) were parameter-specific. These impacts have been summarised in Table 4, which combines conclusions from both the PRCC and Sobol indices.

The RUE was an important driver of the variation in TDW and QS<sub>mobile,GL</sub>. It was the most impacting parameter for TDW with a positive trend that was stable across all S supply conditions (amounts and timing), meaning that increasing its value leads to higher biomass at the onset of pod formation. There is a direct correlation between RUE and TDW from Monteith's model [38], which is used in the SuMoToRI model [2]. In contrast, the two other parameters, β and SLA, had very little impact on TDW, which made sense because they directly relate to the leaves. According to the total Sobol indices, the allocation of carbohydrates to the leaves (β) was slightly more influential than the SLA (with values ranges of 0.072–0.12 for β and 0.041 to 0.094 for SLA). Because no interactions were observed, the effective impact of a given parameter's variation was not influenced by the other parameters, which confirms that they each described distinct processes (i.e. biomass production vs. biomass allocation).

The SAs also revealed that the most impacting parameters on LAI<sub>GL</sub> were β, and to a lesser extent, SLA and RUE. As expected, increased C-allocation to the leaves (β) and SLA favoured greater leaf area expansion while higher RUE had little impact. While the effect of β was stable with increasing S-amounts when S was provided at bolting alone or at bolting and flowering, the effects of SLA and RUE were influenced by the S amounts in these timing conditions. This was consistent with prior observations during the model calibration and evaluation steps that (i) led to different RUE values according to the S-supply conditions [2] and (ii) pointed out that the SLA decreased under S-limiting conditions, meaning that the leaf thickness was greater (Brunel-Muguet, unpublished). This leaf plasticity was showed to increase the photosynthetic activity, thus setting a compensatory mechanism against drastic S limitation, which is known to impair photosynthesis and C metabolism [39].

Regarding the amount of mobile S within the leaves (QS<sub>mobile,GL</sub>), the respective impact of the parameters varied according to S supply conditions. The Sobol indices indicated that RUE was the main driver followed by β and SLA. Increasing RUE, β and SLA led to lower QS<sub>mobile,GL</sub>, which was amplified by increases in S fertilisation when supplied either once at bolting or
twice under the fractioned condition. In contrast, the impact of the three parameters was much lower under all levels of S supply at flowering. As observed, increasing the RUE favoured biomass production and consequently leaf expansion. This in turn increased the leaf S-struc- tural requirements and thus depleted the mobile S amounts stored within the leaves. Another reason that could account for these effects is the dilution of S within the plant. As the size of the source organs increases, the concentration of S is diluted, leading to depletion of S mobile reserves. This effect of increased RUE was observed for optimal timing i.e. at bolting and in the fractioned condition. It was also lower as the S supplies decreased.

## A preliminary approach to design ideotypes adapted to specific S-fertilising strategies

In our study, we made the assumption that ideal parameter combinations combined with ade- quate S-fertilisation strategies should offer the most efficient plant status for achieving the reproductive phase and then high yield. Our analyses revealed that the optimised combina- tions for the best performance at the onset of pod formation (i.e. high biomass and green leaf area to maximise plant development and light interception and high mobile S to allow ade- quate S storage for growing pods) were specific to each output, thus meaning that there will be trade-offs to attain them in combination (Table 3). For instance, the optimal RUE value for the highest TDW (for the RC<sub>opt</sub>) is 4.58, whereas it is 2.59 to attain the highest LAI<sub>GL</sub>. Moreover, in some cases, it is crucial to consider several parameters together because of their strong inter- actions, as illustrated for QS<sub>mobile,GL</sub> (Fig 5B).

Table 4 synthesises the parameter rankings for a given output (i.e. plant biomass, leaf expansion, S storage) according to the S supply strategies. It indicated that the most beneficial parameter for targeting would differ depending on whether we intended to improve plant bio- mass, leaf expansion or S storage within the leaves with respect to the S-fertilisation designs. Indeed, RUE was the driver of plant biomass and its impact was stable regardless of amount and timing of S-fertilisation (stage and fractioning). Leaf expansion was mainly boosted by increased C allocation to the leaves (β) and to a lesser extent by increased specific leaf area and RUE, in all S-fertilisation designs apart from low S amounts (20U) when supplied at bolting (GS30) or in the fractioned condition (GS360+690). Impacts of increased RUE and SLA dif- fered depending on the S-fertilisation designs. While the impact of SLA was higher with increasing rates of S-fertilisers at flowering (GS60), the impact of RUE was more pronounced in the fractioned condition with extreme amounts (20 and 300U).

These SAs revealed that the association between high biomass or leaf area, and S storage for growing was not trivial because the plant parameter drivers were distinct and the intensity of their impact was also dependent on the S fertilisation designs. Genotypes with high radiation efficiency, large light interception surfaces (per biomass unit) or high leaf biomass allocation would benefit from early and fractioned S-fertilisation strategies.

## Towards new S management strategies

These results confirm the need to think about new designs or cultural practices for S fertilisa- tion. First, they highlighted the negative impact of delayed fertilisation (i.e. at flowering). When applied at this stage (GS60), a higher S supply could not compensate for the earlier reductions in vegetative growth. In addition, our results indicated the relevance of fractioning S supply compared to the conventional single supply at bolting, and thus broadened fertilisa- tion schemes. In our conditions, S-fractioning led to similar plant performances at the onset of pod formation. This strategy also prevented excessive sulfate being stored in the leaves and then its potential loss when leaves detached. From an ecological and economic perspective,
fractioning would allow better nutrient adjustment rather than the conventional single supply and thus help in preventing over fertilisation. Nevertheless, our results showed that even under the most optimised plant parameter combination, the increase in TDW of around 13% under the fractioned condition was obtained with the highest S supply, which ultimately moderates the sustainable effect of fractioning. In addition, it is known that S nutrition is closely related to N nutrition in many species like oilseed rape [21,22,40], wheat [41,42] and maize [43]. Fismes et al. [44] reported that in oilseed rape, the balance between the S and N rate determines their use efficiencies, which are synergistic at optimum rates and antagonistic at excessive levels of a single element. S fertilisation is known to improve N use efficiency and to maintain high seeds quality [45], thus highlighting the importance of balanced N:S ratio. In our study, the fifteen S fertilisation conditions were provided with non-limiting N supply. For the lowest S supply conditions, the N:S ratio was not as well balanced as for the conditions with the conventional S supply, which could impact crop performances. N inputs are drivers of this ratio meaning they can modulate crop responses to S supply. In this context, implementing S: N ratio thresholds in the model could help monitor S fertilisation to better monitor effective S assimilation, which is partly determined by N availability.

## Conclusion

A better understanding of the relationship between specific model parameters and crops at the onset of pod formation is essential to adjust cultural practices under contrasting S fertilisation management. For this purpose, we used GSAs to rank driver model parameters including RUE, β and SLA according to their impact on main output variables related to growth, leaf expansion and S remobilisation. In our study, we compared these plant parameters combinations to distinguish S-fertilisation designs. This enabled the best S-fertilisation designs to be determined under specific plant parameters combinations. This approach underscores the advantages of using adapted varieties for a specific fertilisation context. Our study also highlighted the importance of S fertilisation management in terms of the amount supplied or the timing of the application (stage and fractioning). We showed that supplying more than the recommended inputs did not lead to proportional improvements in plant biomass or leaf expansion and that delaying supply until flowering had a strong negative impact. If fractioning designs seem to be competitive with the conventional single supply (75U SO₃ at the end of winter), significant increased performance would require higher S rates. These findings reinforce the need to develop tools to measure plant S status and soil S availability in situ, which would help in adjusting fractioned S supplies throughout the crop cycle.

## Supporting information

**S1 Fig. Climatic data used for model simulation and sensitivity analyses.** Daily mean temperature in °C (A) and daily mean Photosynthetically Active Radiation (PAR) in MJ.m⁻² (B) in Saint-Martin-de-Hinx averaged over 2005 to 2015 (source: CLIMATIK https://intranet.inra.fr/climatik_v2) for the simulated periods. (TIF)

## Acknowledgments

The authors thank Laurence Cantrill for his assistance in editing, Christophe Nguyen for his contribution to the Sobol script and the MEXICO network. We also thank two anonymous reviewers for their comments on the manuscript. This work was performed, in partnership with the SAS PIVERT, within the framework of the French Institute for Energy Transition
(Institut pour la Transition Energétique (ITE) P.I.V.E.R.T. (www.institut-pivert.com) selected as an Investment for the Future (“Investissements d'Avenir”). This work was supported, as part of the Investments for the Future by the French Government under the reference, ANR-001-01.

## Author Contributions

**Conceptualization:** Emilie Poisson, Sophie Brunel-Muguet, François Kauffmann, Jacques Trouverie, Jean-Christophe Avice, Alain Mollier.

**Data curation:** Emilie Poisson, François Kauffmann.

**Formal analysis:** Emilie Poisson, Sophie Brunel-Muguet, François Kauffmann, Alain Mollier.

**Funding acquisition:** Jean-Christophe Avice.

**Investigation:** Emilie Poisson, Sophie Brunel-Muguet, Jacques Trouverie, Jean-Christophe Avice, Alain Mollier.

**Methodology:** Emilie Poisson, Sophie Brunel-Muguet, François Kauffmann, Alain Mollier.

**Project administration:** Emilie Poisson.

**Resources:** Emilie Poisson.

**Software:** Emilie Poisson, François Kauffmann, Alain Mollier.

**Supervision:** Emilie Poisson, Sophie Brunel-Muguet, Jacques Trouverie, Jean-Christophe Avice, Alain Mollier.

**Validation:** Emilie Poisson, Sophie Brunel-Muguet, François Kauffmann, Jacques Trouverie, Jean-Christophe Avice, Alain Mollier.

**Visualization:** Emilie Poisson, Alain Mollier.

**Writing – original draft:** Emilie Poisson, Sophie Brunel-Muguet.

**Writing – review & editing:** Emilie Poisson, Sophie Brunel-Muguet, François Kauffmann, Jacques Trouverie, Jean-Christophe Avice, Alain Mollier.

## References

1. Brisson N, Mary B, Ripoche D, Jeuffroy MH, Ruget F, Nicoulaud B, et al. STICS: a generic model for the simulation of crops and their water and nitrogen balances. I. Theory and parameterization applied to wheat and corn. Agronomie 1998; 18:311–46.

2. Brunel-Muguet S, Mollier A, Kauffmann F, Avice J-C, Goudier D, Sénécal E, et al. SuMoToRI, an Ecophysiological Model to Predict Growth and Sulfur Allocation and Partitioning in Oilseed Rape (Brassica napus L.) Until the Onset of Pod Formation. Frontiers in plant science 2015; 6:993. https://doi.org/10.3389/fpls.2015.00993 PMID: 26635825

3. Jamieson PD, Semenov MA, Brooking IR, Francis GS. Sirius: a mechanistic model of wheat response to environmental variation. European Journal of Agronomy 1998; 8:161–79.

4. Jeuffroy MH, Recous S. Azodyn: A simple model simulating the date of nitrogen deficiency for decision support in wheat fertilization. European Journal of Agronomy 1999; 10:129–44.

5. Mollier A, De Willigen P, Heinen M, Morel C, Schneider A, Pellerin S. A two-dimensional simulation model of phosphorus uptake including crop growth and P-response. Ecological Modelling 2008; 210:453–64.

6. Robertson MJ, Lilley JM. Simulation of growth, development and yield of canola (Brassica napus) in APSIM. Crop and Pasture Science 2016;67.

7. Asseng S, Zhu Y, Wang E, Zhang W. Crop modeling for climate change impact and adaptation. In: Crop Physiology: Applications for Genetic Improvement and Agronomy: Second Edition. Elsevier; 2014:505–46.
8. Challinor AJ, Watson J, Lobell DB, Howden SM, Smith DR, Chhetri N. A meta-analysis of crop yield under climate change and adaptation. Nature Climate Change 2014; 4:287-91.

9. Chenu K, Porter JR, Martre P, Basso B, Chapman SC, Ewert F, et al. Contribution of Crop Models to Adaptation in Wheat. Trends in Plant Science 2017; 22:472-90. https://doi.org/10.1016/j.tplants.2017.02.003 PMID: 28389147

10. Kumar M. Impact of climate change on crop yield and role of model for achieving food security. Environmental Monitoring and Assessment 2016; 188:465. https://doi.org/10.1007/s10661-016-5472-3 PMID: 27418072

11. Nuttall JG, O'Leary GJ, Panozzo JF, Walker CK, Barlow KM, Fitzgerald GJ. Models of grain quality in wheat—A review. Field Crops Research 2017; 202:136-45.

12. Cariboni J, Gatelli D, Liska R, Saltelli A. The role of sensitivity analysis in ecological modelling. Ecological Modelling 2007; 203:167-82.

13. Iooss B, Lemaître P. A Review on Global Sensitivity Analysis Methods. Uncertainty management in Simulation-Optimization of Complex Systems: Algorithms and Applications 2015;101-22.

14. Saltelli a., Tarantola S, Campolongo F. Sensitivity analysis as an ingredient of modeling. Statistical Science 2000; 15:377-95.

15. Colbach N, Molinari N, Clermont-Dauphin C. Sensitivity analyses for a model simulating demography and genotype evolutions with time. Ecological Modelling 2004; 179:91-113.

16. Oenema O, Postma R. Managing Sulphur in Agroecosystems. In: Sulphur in Plants. Dordrecht: Springer Netherlands; 2003. page 45-70.

17. Schnug E, Haneklaus S, Murphy D. Impact of sulphur fertilization on fertilizer nitrogen efficiency. Sulphur in Agriculture 1993; 17:8-12.

18. Brunel-Muguet S, D'Hooghe P, Bataillé M-P, Larré C, Kim T-H, Trouverie J, et al. Heat stress during seed filling interferes with sulfur restriction on grain composition and seed germination in oilseed rape (Brassica napus L.). Frontiers in plant science 2015; 6:213. https://doi.org/10.3389/fpls.2015.00213 PMID: 25914702

19. D'Hooghe P, Escamez S, Trouverie J, Avice JC. Sulphur limitation provokes physiological and leaf proteome changes in oilseed rape that lead to perturbation of sulphur, carbon and oxidative metabolisms. BMC Plant Biology 2013; 13:23. https://doi.org/10.1186/1471-2229-13-23 PMID: 23391283

20. Dubousset L, Etienne P, Avice JC. Is the remobilization of S and N reserves for seed filling of winter oilseed rape modulated by sulphate restrictions occurring at different growth stages? Journal of Experimental Botany 2010; 61:4313-24. https://doi.org/10.1093/jxb/erq233 PMID: 20693411

21. Dubousset L, Abdallah M, Desfeux AS, Etienne P, Meuriot F, Hawkesford MJ, et al. Remobilization of leaf S compounds and senescence in response to restricted sulphate supply during the vegetative stage of oilseed rape are affected by mineral N availability. Journal of Experimental Botany 2009; 60:3239-53. https://doi.org/10.1093/jxb/erp172 PMID: 19553370

22. Abdallah M, Dubousset L, Meuriot F, Etienne P, Avice JC, Ourry A. Effect of mineral sulphur availability on nitrogen and sulphur uptake and remobilization during the vegetative growth of Brassica napus L. Journal of Experimental Botany 2010; 61:2635-46. https://doi.org/10.1093/jxb/erq096 PMID: 20403880

23. Girondé A, Dubousset L, Trouverie J, Etienne P, Avice J-C. The impact of sulfate restriction on seed yield and quality of winter oilseed rape depends on the ability to remobilize sulfate from vegetative tissues to reproductive organs. Frontiers in plant science 2014; 5:695. https://doi.org/10.3389/fpls.2014.00695 PMID: 25566272

24. Morris MD. Factorial Sampling Techniques for Preliminary Computational Experiments. Technimetrics 1991; 33:161-74.

25. Iooss B. Revue sur l'analyse de sensibilité globale de modèles numériques. Journal de la société française de statistique, Société française de statistique et société mathématique de France 2011; 152:1-23.

26. Jullien A, Allirand JM, Mathieu A, Andrieu B, Ney B. Variations in leaf mass per area according to N nutrition, plant age, and leaf position reflect ontogenetic plasticity in winter oilseed rape (Brassica napus L.). Field Crops Research 2009; 114:188-97.

27. Randall PJ, Wang Q, Hocking PJ, Pinkerton A. Critical values for sulfur in young plants of oilseed rape (Brassica napus L.) determined with reference to dry weight, leaf area and specific leaf weight. In: Plant Nutrition for Sustainable Food Production and Environment. Dordrecht: Springer Netherlands; 1997. page 335-9.

28. Khan NA, Mobin M, Samiullah. The influence of gibberellic acid and sulfur fertilization rate on growth and S-use efficiency of mustard (Brassica juncea). Plant and Soil 2005; 270:269-74.
29. Poisson E, Mollier A, Trouverie J, Avice J-C, Brunel-Muguet S. SuMoToRI model simulations for optimizing sulphur fertilization in oilseed rape in the context of increased spring temperatures. European Journal of Agronomy 2018; 97:28-37.

30. Ingestad T, Lund A-B. Theory and techniques for steady state mineral nutrition and growth of plants. Scandinavian Journal of Forest Research 1986; 1:439-53.

31. Oscarson P, Ingemarsson B, Larsson C-M. Growth and nitrate uptake properties of plants grown at different relative rates of nitrogen supply. II. Activity and affinity of the nitrate uptake system in Pisum and Lemna in relation to nitrogen availability and nitrogen demand. Plant, Cell and Environment 1989; 12:787-94.

32. Morrison MJ, Mcvettyi PBE, Shaykewich CF3. The determination and verification of a baseline temperature for the growth of westar summer rape. Can. J. Plant Sci. Can. J. Plant Sci 1989; 69:455-64.

33. Lancashire PD, Bleiholder H, Van Den Boom T, Langelüddeke P, Stauss R, Weber E, et al. A uniform decimal code for growth stages of crops and weeds. Annals of Applied Biology 1991; 119:561-601.

34. Mckay MD, Beckman RJ, Conover WJ. Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code. Technometrics 1979; 21:239-45.

35. Saltelli A. Sensitivity analysis for importance assessment. In: Risk Analysis. Blackwell Publishing, Inc.; 2002. page 579-90. PMID: 12088235

36. Sobol IM. Sensitivity analysis for nonlinear mathematical models. Matematicheskoe Modelirovanie 1990; 2:112-8.

37. Team R core. The R Project for Statistical Computing. http://www.R-Project.Org/ 2013;1-12.

38. Monteith JL, Moss CJ. Climate and the Efficiency of Crop Production in Britain [and Discussion]. Philosophical Transactions of the Royal Society B: Biological Sciences 1977; 281:277-94.

39. D'Hooghe P, Dubousset L, Gallardo K, Kopriva S, Avice J-C, Trouverie J. Evidence for proteomic and metabolic adaptations associated to alterations of seed yield and quality in sulphur-limited Brassica napus L. Molecular & cellular proteomics: MCP 2014; 13:1165-83. https://doi.org/10.1074/mcp.M113.034215 PMID: 24554741

40. Ahmad G, Jan A, Arif M, Jan MT, Khattak RA. Influence of nitrogen and sulfur fertilization on quality of canola (Brassica napus L.) under rainfed conditions. Journal of Zhejiang University SCIENCE B 2007; 8:731-7. https://doi.org/10.1631/jzus.2007.B0731 PMID: 17910116

41. Salvagiotti F, Castellarín JM, Miralles DJ, Pedrol HM. Sulfur fertilization improves nitrogen use efficiency in wheat by increasing nitrogen uptake. Field Crops Research 2009; 113:170-7.

42. Salvagiotti F, Miralles DJ. Radiation interception, biomass production and grain yield as affected by the interaction of nitrogen and sulfur fertilization in wheat. European Journal of Agronomy 2008; 28:282-90.

43. Dawar K, Iqbal J, Wahab S. Effect of Sulfur on Nitrogen Use Efficiency and Yield of Maize Crop. Advances in Environmental Biology 2016; 10:85-90.

44. Fismes J, Vong PC, Guckert A, Frossard E. Influence of sulfur on apparent N-use efficiency, yield and quality of oilseed rape (Brassica napus L.) grown on a calcareous soil. European Journal of Agronomy 2000; 12:127-41.

45. Anjum NA, Gill SS, Umar S, Ahmad I, Duarte AC, Pereira E. Improving growth and productivity of Oleiferous Brassicas under changing environment: significance of nitrogen and sulphur nutrition, and underlying mechanisms. TheScientificWorldJournal 2012; 2012:657808. https://doi.org/10.1100/2012/657808 PMID: 22629181