Analysis of the Similarity between in Silico Ideotypes and Phenotypic Profiles to Support Cultivar Recommendation—A Case Study on Phaseolus vulgaris L.

Abstract

Cultivar recommendation is a key factor in cropping system management. Classical approaches based on comparative multi-environmental trials can hardly explore the agro-climatic and management heterogeneity farmers may have to face. Moreover, they struggle to keep up with the number of genotypes commercially released each year. We propose a new approach based on the integration of in silico ideotyping and functional trait profiling, with the common bean (Phaseoulus vulgaris L.) in Northern Italy as a case study. Statistical distributions for six functional traits (light extinction coefficient, radiation use efficiency, thermal time to first pod and maturity, seed weight, plant height) were derived for 24 bean varieties. The analysis of soil, climate and management in the study area led us to define 21 homogeneous contexts, for which ideotypes were identified using the crop model STICS (Simulateur mulTIdisciplinaire pour les Cultures Standard), the E-FAST (Extended Fourier Amplitude Sensitivity Test) sensitivity analysis method, and the distributions of functional traits. For each context, the 24 cultivars were ranked according to the similarity (weighted Euclidean distance) with the ideotype. Context-specific ideotypes mainly differed for phenological adaptation to specific combinations of climate and management (sowing time) factors, and this reflected in the cultivar recommendation for the different contexts. Feedbacks from bean technicians in the study area confirmed the reliability of the results and, in turn, of the proposed methodology.

1. Introduction

Crop performance is the result of the complex interactions between genotypic factors and the agronomic and environmental context in which the crop is grown. Understanding the genotype (G) by environment (E) and by management (M) interactions (G × E × M, [1]) is thus crucial for cultivar recommendation, as it allows exploring the opportunities or disadvantages of different genotypes in a specific cultivation area [2].

Multi-environmental trials (METs), where the performances of promising genotypes are evaluated in different environments [3], represent one of the main tools to account for phenotypic plasticity within cultivar recommendations [4]. In fact, the analysis of crop responses across the target population of environments (TPE, i.e., the set of environmental and agronomic conditions in which the varieties produced by a breeding program will be grown [5]) allows us to identify the best genotypes for a given agro-environmental context [6].

Although METs can provide sound information about the genotype potential in the TPEs, their effectiveness is often limited by the amount of data needed to actually capture the effects of the variability in the environmental and management factors of the target area [7,8], which may raise the cost of field trials to unsustainable levels. As a result, METs can be biased by non-representative sampling of the TPEs [9], which may decrease the predictive capability for environments and cropping systems not—or insufficiently—represented in the field testing phase [10].

Given their capability to reproduce genotype growth and development as a response to environmental and management factors [11], a variety of modelling approaches have been developed to better understand G × E × M interactions in multi-site field trials and to increase the effectiveness of the way environments are sampled within METs. Crop models have been used to broadly characterize cropping systems in terms of water-deficit patterns, by identifying drought intensity, seasonal dynamics and frequency of occurrence over long time series and large areas [9,12,13], ultimately to support germplasm testing strategies. Model-based approaches to directly assist variety evaluations have also been proposed [4], which use data collected in national-scale METs for cultivar-specific model calibration. The model was then used to reproduce phenotypic plasticity over the entire target area, thus amplifying the range of agro-climatic conditions explored while testing the varieties. The use of crop models to analyze G × E × M interactions also expanded towards the definition of context-specific ideotypes [14,15], which can be defined as an ideal combination of functional traits [16] expected to perform in a desired manner within a given agro-climatic context. Plant morphological and physiological traits can be codified as model parameters, which can be used to evaluate the impact of varying one or more plant traits on crop growth and development in a wide range of soil, climate, and management conditions, even including entire production areas [17].

In those terms, model-based ideotyping and cultivar recommendation share the goal of identifying a plant type able to reach specific productive goals in a given agro-climatic context, leading to the hypothesis that the use of in-silico ideotyping techniques could allow a more comprehensive evaluation of G × E × M interactions while supporting cultivar recommendations [18].

To test this hypothesis, we focused on the common bean (Phaseolus vulgaris L.) as a case study, for its importance as a major source of dietary protein in both less developed and developing countries, and for its soil-enriching role in cropping systems [19]. Among the major food crops, beans show one of the highest levels of variation in growth habits, seed characteristics and adaptation capabilities [20,21], which turn into a huge number of accessions available worldwide (almost 40,000 in the sole International Centre for Tropical Agriculture—CIAT—gene bank https://www.genebanks.org/resources/crops/beans/). From a long-term perspective, this germplasm diversity represents a precious resource for genetic improvement targeting higher yields and resistance/tolerance to biotic and abiotic stressors [22]. In the short term, this makes cultivar choice a key step to take full advantage of the large variability in cultivar features, ultimately to obtain high and stable yields in any given cultivation environment [23].

The objectives of this study were (i) developing a new methodology for supporting cultivar recommendation based on the similarity between in silico ideotypes and functional trait profiles of available varieties and (ii) testing the methodology on a real case study. The approach developed is generic, and is thus easily transferable to other crops and traits.

The crop model STICS (Simulateur mulTIdisciplinaire pour les Cultures Standard) [24,25] was used as a simulation engine, this being—according to the authors’ knowledge—the first time this model has been tested for the common bean (Phaseouls vulgaris L.).

2. Materials and Methods

2.1. Characterization of the Study Area

The Emilia-Romagna region (Northern Italy) was selected as a study area, given its primary role in terms of national bean production and harvested areas. The region is characterized by intensive agricultural production systems and by a marked variability in terms of environmental conditions, due to the presence of coastal zones in the eastern part of the region and of the Appennini mountains in the south-west side (Figure 1).

Different contexts were identified in the study area, according to the variability in soil, climate and management conditions, this being the first step to derive a reliable system to support cultivar recommendations [26]. For the zoning process, only the areas with an altitude lower than 500 m a.s.l. were considered, given that horticulture is performed only in the flat areas of the region.

The soil texture map of Emilia-Romagna was retrieved from the regional geological service (https://geo.regione.emilia-romagna.it/geocatalogo/), and three soil classes (fine: clay content >35% and sand content <65%; coarse: clay content <17.5% and sand content >65%; medium: all remaining soils [27]) were adopted to derive zones homogeneous in terms of soil properties (Figure 2a).

Variability in management factors was accounted for by considering two sowing dates, representative of the main sowing periods adopted by farmers for dry bean production: 15th April for early sowings and 15th May for late sowings (early and late sowings hereafter). No differences across the region in terms of water and nitrogen management were considered, the crop being always fully irrigated and with N requirements satisfied by about 50 kg N ha⁻¹ provided as urea (standard practices in the area) in addition to the N made available by symbiotic fixation.

Daily weather data from 2009 to 2018 were retrieved from the high-resolution (5 km × 5 km) gridded climate dataset available at the Regional Agency for Prevention, Environment and Energy of the Emilia-Romagna website (ARPAE), which is based on 60 weather stations for temperature and 254 for rainfall [28]. Overall, the study area consists of 679 cells of 5 km × 5 km. For each of them, global solar radiation (RAD, MJ m⁻² d⁻¹) and reference evapotranspiration (ET0, mm d⁻¹) were derived by using the approaches reported by [29,30].

For each combination of climate cell × sowing date, the Synthetic AgroMeteorological indicator (SAM; −1 to +1; Equation (1); [31]) was calculated as:

$$SAM=\frac{Rain-ET0}{Rain+ET0}$$

(1)

where Rain is the 10-year average cumulated rainfall (mm) in the 100 days after sowing (average duration of the cropping season for bean in Northern Italy) and ET0 is the ET0 calculated over the same period.

Four classes were used to cluster the variability observed in the SAM indicator, leading to the identification of three and four climate subzones for early and late sowings, respectively (Figure 2b).

The intersection of the climate and soil zonations led to the identification of 21 homogeneous contexts (Figure 2c; Table S1), resulting from nine combinations of soil and climate classes for early sowings, and twelve combinations for the late ones.

2.2. Calibration of the Crop Model

The crop model STICS [24,25,32] was selected for the study, given its proven reliability for the simulation of legume species [33,34] and its high level of plasticity [35]. Model plasticity, which can be defined as “the model aptitude to change the sensitivity to parameters while changing the conditions explored” [36], is indeed a desirable model feature when G × E × M interactions need to be explored.

STICS is a generic soil-crop model, reproducing crop growth and development with a daily time step and simultaneously estimating the carbon (C), nitrogen (N), and water balances in the soil–crop system. It accounts for both determinate and indeterminate growth habits, the latter reproduced by considering the trophic interactions between multiple cohorts of fruits. Biomass accumulation is reproduced with a net photosynthesis approach based on radiation use efficiency, as modulated by the effect of temperature, water stress, N availability, excess of radiation, and concentration of CO₂ in the atmosphere. Canopy light harvesting depends on leaf area index (LAI) and on the light extinction coefficient (Lambert-Beer’s approach), with fruit pods also contributing to light interception. Crop development is reproduced based on thermal time accumulation, computed as a function of daily mean crop temperature and parameters representing base, optimum, and critical temperatures, and modulated by the effect of vernalization requirements, photoperiod response and water stress. LAI dynamics are simulated by considering development stage and the crop response to temperature, plant density, water and nitrogen stress. A complete description of the model algorithms is available in the seminal literature and in the model documentation [24,25].

In order to adapt the model to the bean varieties used in the study area, parameters were calibrated using data collected in dedicated field experiments carried out between 2016 and 2018 at five sites. The two common bean cultivars Etna and Taylor’s horticultural (hereafter referred to as Taylor’s) were grown, given they are among the most widely cultivated in the study area. Etna was sown on 11 April 2017 in Guardamiglio (45.13° N, 9.64° E, 49 m a.s.l.), on 20 April 2017 in Castelvetro Piacentino (45.08° N, 9.96° E, 39 m a.s.l.), and on 8 May 2018 in Cadriano (44.55° N, 11.41° E, 32 m a.s.l.), whereas Taylor’s was sown on 26 April 2016 in Jolanda di Savoia (44.88° N, 11.97° E, −1 m a.s.l.), on 12 April 2017 in Gossolengo (45.01° N, 9.62° E, 86 m a.s.l.), and on 8 May 2018 in Cadriano. Soil texture was clay-loam in Guardamiglio and Castelvetro Piacentino, clay in Jolanda di Savoia, and silty-clay-loam in Cadriano and Gossolengo (USDA classification), with soil organic carbon ranging from 1.3% (Cadriano) to 5% (Jolanda di Savoia). Field management was allowed to keep the crop weed-, disease- and pest-free, with an optimal water and nutrient supply.

In the three field experiments carried out in 2017, phenology was recorded by using the BBCH scale for bean [37], and dry aboveground biomass (AGB, t ha⁻¹), leaf area index (LAI, -), specific leaf area (SLA, m² kg⁻¹), canopy height (cm), and the number of leaves and pods per plant were measured at three random points in the field in four sampling events, i.e., at the beginning of flowering (BBCH 60), the beginning of pod development (BBCH 71), the beginning of ripening (BBCH 81), and at maturity (BBCH 89). For AGB, twenty plants for each point were sampled and dried until a constant weight was found. Dry biomass of pods and, at maturity, of seeds was also measured. LAI was estimated by using the PocketLAI smart-app [38], with five measurement replicates at each point. SLA was derived by digitalizing sampled leaves from three random plants and calculating the leaf area to dry leaf mass ratio. In the 2016 field experiment (Jolanda di Savoia), the same sampling scheme as described above was adopted but LAI was not measured. In the experiment in Cadriano, LAI and AGB were measured for both Etna and Taylor’s in two moments, at inflorescence emergence (BBCH 51) and at the start of pod development (BBCH 69). At maturity, dry grain yield was measured.

The STICS parametrization for bean was developed by using a trial-and-error approach, targeting the highest agreement between observed and simulated values for phenology, LAI, aboveground biomass, pod biomass, and yield. For both calibration and validation, model performance was evaluated by using five indices of agreement: mean absolute error (MAE [39]; from 0 to + ∞, optimum: 0); relative root mean square error (RRMSE [40]; from 0 to + ∞, optimum: 0); Nash and Sutcliffe modelling efficiency (EF [41]; from −∞ to 1, optimum: 1); coefficient of residual Mass (CRM [42]; −∞ to +∞, optimum: 0; positive in case of underestimation and vice versa); and coefficient of determination (R²; from 0 to 1, optimum: 1) of the linear regression between observed and simulated data.

2.3. Determination of Functional Trait Values

Twenty-four bean genotypes—including 5 commercial cultivars and 19 new varieties (Table S2)—were evaluated in a dedicated field experiment carried out in Cadriano (44.55° N, 11.41° E) in 2018. Plants were sown on 8 May in 3 m × 7 m plots with a density of 30 plants m⁻² and grown under optimal conditions (no limiting factors other than radiation and temperature). The genotypes—all with a determinate growth habit and not sensitive to photoperiod—can be regarded as a representative of the borlotto beans cultivated in Italy, since they cover most of the harvested area for this crop.

Six functional traits were considered related to phenological development (thermal time to begin pod development and thermal time to reach maturity), canopy structure (light extinction coefficient and plant height), photosynthetic efficiency (radiation use efficiency), and yield formation (maximum seed weight).

Dates of phenological stages were recorded at complete emergence (BBCH 10), beginning of pod development (BBCH 69), and full ripening (BBCH 89), and they were then used to derive the thermal time requirements to begin pod development and maturity, using 7 °C as base temperature.

The light extinction coefficient was estimated at first pod appearance by using the smart app PocketPlant3D [43]. PocketPlant3D allows obtaining a 3D distribution of the angles of photosynthetic tissues, which is used to estimate the light extinction coefficient according to Campbell [44]. The 3D scanning was performed on five plants per plot. At the same time, crop height was measured by averaging five measurements per plot.

Radiation use efficiency (g MJ⁻¹, Equation (2)) was estimated according to Paleari et al. [45] in the time interval between complete canopy closure and the beginning of pod development, in order to avoid uncertainty in RUE estimation due to either incomplete soil coverage or leaf senescence:

$$RUE=\frac{AGB}{0.5·RAD·\left(1-e^{\left(-k·LAI\right)}\right)}·\frac{1}{Tlim}$$

(2)

where AGB is the total aboveground biomass increase (g m⁻²) in the time interval Δt estimated from the dry biomass measured on five plants per plot at BBCH 51 (flower buds visible) and BBCH 69 (first pod appearance); 0.5 is a conversion factor to derive photosynthetic active radiation (PAR) from global solar radiation; RAD is the cumulated global solar radiation (MJ m^−²) during Δt; k is the extinction coefficient for solar radiation estimated for each cultivar; LAI is the mean leaf area index measured during the time interval Δt using the smart app PocketLAI [38]; Tlim is the mean thermal limitation to net photosynthesis during Δt, which was estimated as a function of cardinal temperatures (trapezoid response curve parameterized with 7 °C, 28 °C, 32 °C, and 42 °C) and mean daily air temperature [25].

At maturity, 50 plants per plot were sampled for yield determination, and seed weight was derived as the average weight of a random sample of 100 dry seeds.

The Grubbs [46] test was applied to detect the presence of outliers for all measurements.

Once the values of the six functional traits for the 24 varieties were determined, their normality was checked using the Shapiro–Wilk test [47], and distribution parameters (mean and standard deviation) were derived for each trait.

2.4. Ideotype Design and Cultivar Recommendation

The distributions of trait values were used to derive ideotypes specific for each of the 21 simulation contexts (homogeneous for soil, climate, and management) using the methodology proposed by Paleari et al. [14]. Under the assumption of a close relationship between model parameters and plant traits (

Table 1. Traits considered in the study, with corresponding STICS parameters and statistical distribution (mean ± standard deviation) derived from the characterization of 24 bean varieties. Results of the Shapiro–Wilk normality test are also shown.

Process	Trait	Parameter (Acronym)	Unit	Distribution	Shapiro-Wilk Test
Photosynthesis	Radiation use efficiency	Maximum radiation use efficiency of the vegetative phase (efcroiveg)	g MJ−1	Normal (2.2 ± 0.5)	W = 0.949, p-value = 0.256
Yield formation	Seed weight	Maximum grain weight at 0% humidity (pgrainmaxi)	g	Normal (0.46 ± 0.04)	W = 0.987, p-value = 0.983
Canopy structure	Plant height	Maximum plant height (hautmax)	m	Normal (0.6 ± 0.09)	W = 0.921, p-value = 0.062
	Light extinction coefficient	Extinction coefficient for solar radiation (extin)	-	Normal (0.9 ± 0.05)	W = 0.912, p-value = 0.039 a
Phenological development	Thermal time to first pod	Cumulative thermal time from emergence to the onset of filling of harvested organs (stlevdrp)	°C-day	Normal (530 ± 52)	W = 0.955, p-value = 0.349
	Thermal time to maturity	Cumulative thermal time from beginning of harvested organs filling to maturity (stdrpmat)	°C-day	Normal (761 ± 80)	W = 0.973, p-value = 0.752

^a: slightly deviating from normality but considered as normal for sensitivity analysis given it deviated more from other distributions.

), this approach makes use of crop modelling and sensitivity analysis (SA) techniques to identify the combination of functional trait values (i.e., the ideotype) that maximizes the value of a given objective function (e.g., yield or resources use efficiency).

In this study, the variance-based global SA method Extended Fourier Amplitude Sensitivity Test (E-FAST [48]) was used to explore the parameter hyperspace. The number of combinations was set to 1494 according to Confalonieri et al. [49] This method allows quantifying the fraction of model output variance due to changes in the values of input factors—in this case model parameters—in terms of first and total order effects, the latter including the effect of interactions among parameters. The output variable used for the computation of SA metrics is the composite index (Y_index, Equation (3)) proposed by Ravasi et al. [35], which allows to account for both yield and yield stability:

$$Y_{index}=\left[\left(\frac{Y_i}{Y_{max}}\right)·0.7\right]+\left[\left(1-\frac{C{V}_i}{C{V}_{max}}\right)·0.3\right]$$

(3)

where Y_i and CV_i are the mean and the coefficient of variations of the yield values simulated with the combination of parameters i over the considered time frame (from 2009 to 2018); Y_max and CV_max are the mean and the coefficient of variations of the yield values simulated with the combination of parameters that achieved the highest yield. Y_index was used to rank the combinations of parameters sampled using the E-FAST method and to define the ideotype trait values as the mean of parameter values of the best 1% combinations [35].

For each of the 21 homogeneous contexts, the similarity between the ideotype and the 24 varieties for which functional trait values were determined was quantified using the weighted Euclidean distance method proposed by Carvalho et al. [50] for index-based selection.

This approach is based on the definition of the optimum range, the optimum phenotypic value and the relative importance of each trait as derived according to the breeder’s experience. The penalization of trait values outside the optimum range is achieved by considering that:

$$\mathrm{If}L{I}_j\le X_{ij}\le L{S}_j\mathrm{then}{Y}_{ij}=X_{ij}$$

(4)

$$\mathrm{If}{X}_{ij}<L{I}_j\mathrm{then}{Y}_{ij}=X_{ij}+V{O}_j-L{S}_j;$$

(5)

$$\mathrm{If}{X}_{ij}>L{S}_j\mathrm{then}{Y}_{ij}=X_{ij}+V{O}_j-L{I}_j;$$

(6)

where X_ij and Y_ij are, respectively, the measured and the transformed phenotypic value of the i_th variety for the j_th trait; LI_j, LS_j, and VO_j are, respectively, the lower and the upper limit of the optimum range and the optimum value of the j_th trait.

In this study, we defined the optimum range of each trait as the range of the corresponding parameter in the best 1% of combinations used to derive the ideotype, and the optimum value of the trait (VO_j) as that of the ideotype defined (see Table 1 for the trait-parameter correspondence). Given ideotypes were identified for each of the 21 agro-climatic contexts, LI_j, LS_j, and VO_j were defined for each context too.

The transformed and the optimum phenotypic values of the i_th variety for the j_th trait (Y_ij and VO_j) were then standardized and weighted as:

$$y_{ij}=\sqrt{a_j}·\frac{Y_{ij}}{s\left(Y_j\right)}$$

(7)

$$v{o}_j=\sqrt{a_j}·\frac{V{O}_j}{s\left(Y_j\right)}$$

(8)

where s(Y_j) is the standard deviation of the transformed phenotypic value (Y_j) of all the varieties including the ideotype, and a_j is the relative importance of the j_th trait.

In this study, the trait value of the ideotype used to calculate s(Y_j) was that of the corresponding model parameter (Table 1), and its sensitivity index (total order effect) was used as the value of a_j. This allowed to assign a value of the weighting factor aj that is specific for each of the 21 contexts, thus taking into account the fact that the relative importance of each trait may vary according to soil, climate and management factors.

The index used to score each variety (IDI_i) was then calculated according to the following equation (Equation (9)):

$$ID{I}_i=\sqrt{\frac{1}{n}\times {\displaystyle \sum }_{j=1}^n{\left(y_{ij}-v{o}_j\right)}^2}$$

(9)

The lower the IDI_i, the higher the similarity between a variety and the ideotype. Therefore, for each of the 21 agro-climatic contexts, the variety with the highest similarity with the ideotype was considered to be the recommended one.

3. Results

3.1. Model Evaluation

The calibration of model parameters allowed to successfully reproduce bean growth and development under the conditions explored in the study area, and a satisfying performance was confirmed for the validation datasets (

Table 2. Agreement between observed and simulated data for the calibration and validation phase. MAE (same units of the variable): mean absolute error (from 0 to +∞, optimum: 0); RRMSE (%): relative root mean square error (from 0 to + ∞, optimum: 0); EF (−): Nash and Sutcliffe modelling efficiency (from −∞ to 1, optimum: 1); CRM (−): coefficient of residual mass (from −∞ to +∞, optimum: 0); R² (−): coefficient of determination (from 0 to 1, optimum: 1).

Activity	Variable	MAE	RRMSE	EF	CRM	R2
Calibration	Leaf area index (−)	0.70	29.8	0.38	0.04	0.68
	Aboveground biomass (t ha−1)	0.91	21.4	0.93	0.01	0.99
	Pod biomass (t ha−1)	0.73	42.6	0.77	−0.35	0.94
	Yield (t ha−1)	0.40	12.9	0.10	−0.10	0.64
Validation	Leaf area index (−)	0.82	29.1	0.79	0.12	0.96
	Aboveground Biomass (t ha−1)	0.40	13.0	0.96	0.03	0.98
	Pod biomass (t ha−1)	0.38	29.6	0.93	−0.09	0.95
	Yield (t ha−1)	0.32	8.4	0.31	−0.04	0.69

; Figure 3).

Mean RRMSE was equal to 23.4% (including calibration and validation datasets), with values ranging from 10.7% (yield) to 36.1% (pod biomass), and EF was always positive, with the best results obtained for AGB for both the calibration and validation phase (mean value equal to 0.95). CRM revealed a systematic slight model underestimation for LAI and AGB and an overestimation—although negligible—for pod biomass and yield. Concerning R², it was equal to 0.80 in the calibration dataset (average of all variables), with the lowest value achieved for yield (0.64) and the highest for AGB (0.99), whereas it improved in the validation phase (average value of 0.9), with values ranging from 0.69 (yield) to 0.98 (AGB).

3.2. Variability of Functional Trait Values among Varieties

The analysis of the values of functional traits revealed that a normal distribution properly described the variability observed among the 24 varieties (Figure 4, Table 1). A marked variation in the phenotypic response was found for most of the traits, especially for radiation use efficiency (coefficient of variation was equal to 23.4%, values ranging from 1.4 g MJ⁻¹ to 3.1 g MJ⁻¹). For traits involved with thermal time accumulation, the coefficient of variation was around 10% (9.8% and 10.5% for, respectively, thermal time to first pod and to maturity), and the range of values explored was from 414 °C-day to 607 °C-day for thermal time to first pod, and from 617° C-day to 957 °C-day for thermal time to maturity. A similar variability was found for seed weight (coefficient of variation = 8.3%, values ranging from 0.38 g to 0.55 g) and plant height (coefficient of variation = 13.9%, range was 0.5 m to 0.8 m). Light extinction coefficient was the trait with the lowest heterogeneity among the varieties analyzed, with a coefficient of variation of 5% and values ranging from 0.79 to 1.

3.3. Ideotypes Derived for the Agro-Climatic Contexts

Results of the SA showed that thermal time to first pod was the trait with the largest impact on yield and yield stability for all 21 contexts. The corresponding parameter (stlevdrp) was always ranked first (Figure S1) and explained more than 50% of the variability in the composite index used as a target output variable for the SA (Y_index, Equation (3)). Radiation use efficiency during the vegetative phase (parameter efcroiveg) also played a relevant role, being always ranked second (average total order effect equal to 0.4). The third, fourth and fifth place in the rank were instead more variable, with light extinction coefficient, seed weight and thermal time to maturity (parameters hautmax, pgrainmaxi and stdrpmat, respectively) re-ranking according to the context explored. Thermal time to maturity turned out to be particularly relevant for late sowings and dry conditions (SAM class A and B) and less important in other environments, thus showing a marked heterogeneity in the fraction of Y_index variance (values ranging from 2.1% to 12.5%). Light extinction coefficient and seed weight, instead, did not show any particular pattern, explaining on average 4% and 4.9% of variability in Y_index, respectively. The trait with the lowest impact on Y_index was plant height (parameter hautmax), which accounted for less than 3% of the variability of yield and yield stability regardless of the conditions explored.

The primary role of the traits radiation use efficiency and thermal time to first pod reflected in the features of the ideotypes identified for the 21 contexts (Figure S2). All the ideotypes shared indeed an increased photosynthetic efficiency (+32.4% on average; parameter efcroiveg) and a delay in first pod appearance (+12.4% of thermal time needed to reach first pod on average; parameter stlevdrp), although with a larger variability in the parameter relevance in case of late sowings (Figure S2, lower panel). Reduced plant height and a slight increase in light extinction coefficient were the other features shared among the ideotypes identified for the different contexts (−9.4% and +4% on average, respectively, for the parameters hautmax and extin). Differences among ideotypes were mainly due to changes in thermal time to maturity (parameter stdrpmat), for which an increase was suggested under dry conditions (SAM classes A and B) and a slight decrease under wet ones. This pattern was observed for both early and late sowings, although it was clearer for the latter. The larger heterogeneity in ideotype features in case of late sowings was due to the wider range of conditions explored by the crop in terms of both thermal regimes and relationships between rainfall and evapotranspiration.

3.4. Cultivar Recommendation

Figure 5 and Figure 6 show the similarity between the ideotypes identified for each of the 21 contexts and the 24 varieties profiled in this study. The varieties at the top of the charts in Figure 5 and Figure 6 are the most similar to the context-specific ideotypes, thus being those recommended for the different contexts.

The variability observed among the ideotypes defined for each condition (Figure S2) was reflected in the heterogeneity among variety rankings. In case of early sowings (Figure 5), the high level of convergence observed among the ideotypes—especially for the traits with the highest sensitivity index (i.e., radiation use efficiency and thermal time to first pod, Figure S2, upper panel)—turned into the recommendation of the same variety (Etna) in all the agro-climatic contexts. A slight variability was observed from the second place downwards (Figure 5), although the differences in the similarity indexes between the first- and second-ranked varieties (higher than 20%) suggested an avoidance of recommending varieties other than Etna regardless of the context considered.

The varieties to be recommended in case of late sowings were instead different while moving across the agro-climatic contexts (Figure 6). The variability achieved among the ideotypes identified for alternative climate and soil combinations (Figure S2, lower panel) turned indeed into differences in the varieties to be recommended. Etna was still recommended in most cases (top-ranked according to the similarity with the ideotype in about 80% of the cases), but varieties Taylor’s, Meccano, 18B988, 17B916 and Magico played a key role too. Taylor’s was ranked first in two contexts out of twelve (SAM class A with fine and medium soils), and resulted in a good second choice variety (the difference in the similarity index compared to the top-ranked variety was lower than 20%) in 50% of the contexts (blue dots in Figure 6). Varieties 18B988 and Meccano were the most similar to the ideotype in one context out of twelve (18B988 for fine soil and SAM class C, Meccano for coarse soil and SAM class B), and they were in the pool of second choice varieties in 50% (18B988) and 16% (Meccano) of the cases. The placement of 17B916 and Magico in the rankings was instead a clear example of the large variability observed for late sowings. They resulted in nearly being recommendable (17B916 for fine soil and SAM class C, Magico for medium soil and SAM class A) in one context out of twelve, but then they played a marginal role—down to the tenth position in the ranking for Magico—in all the other contexts.

Another difference between the results obtained for early and late sowings was that for the latter it was possible to identify second choice varieties (blue in Figure 6, i.e., genotypes with an intermediate position between recommended cultivar (green in Figure 6) and those regarded as not suitable (red) for a specific growing contexts), thus presenting more options to farmers for variety choice. Differences between early and late sowings were found also in terms of the number of varieties recommended (green in Figure 5 and Figure 6). In the case of late sowings, indeed, up to three varieties were found to be very similar to the ideotype defined for fine soil and SAM class C, whereas only one cultivar resulted in being recommendable for early sowings regardless of the conditions explored (Figure 5).

Both early and late sowings presented a smooth worsening of the similarity index while moving from the top to the bottom of the charts shown in Figure 5 and Figure 6 until the fourth last position in the rank, where a steep decline in the similarity index was observed. The four bottom-ranked varieties were the same for all the 21 contexts analyzed, thus suggesting poor performance for these genotypes regardless of the sowing time, soil properties and climate conditions.

4. Discussion

The distributions derived for the six functional traits highlighted the marked heterogeneity in the pool of bean varieties, although the distribution means were in line with those reported in the literature for the common bean. For instance, the values of radiation use efficiency and light extinction coefficient estimated in this study were consistent with that reported by other authors [51,52] for common beans grown under optimal conditions. In the same way, the mean values for plant height, seed weight, and phenological traits were similar to what was reported for Italian common bean genotypes by Mercati et al. [53] and Campion et al. [54]. Despite distributions being derived using data from a single season, they can be considered as reliable. Indeed, the optimal growing conditions (for water, nutrients, weeds, and pests) and the methodology adopted to account for the effects of temperature on traits affected by thermal limitation (i.e., radiation use efficiency [25]) allowed a reduction of the effect of specific growing conditions on trait value determination [45]. Moreover, the phenotyping activities performed in this study focused on functional traits, which are less affected by G × E × M interactions compared to performance traits [16]. However, the purpose of this study was to present a new methodology to support cultivar recommendations. Therefore, different approaches for the characterization of pools of genotypes can be easily incorporated in the workflow we proposed, e.g., by extending the methodology to traits involved with tolerance/resistance to biotic and abiotic stressors, which can also be accounted for in ideotyping studies (e.g., [14]). The method can be applied also in cases of large variability among genotypes, given that this variability will reflect in the trait distributions used for sensitivity analysis and thus in the variability among the combinations of trait values used to design the ideotypes.

According to the author’s knowledge, this is the first time the crop model STICS is used for the common bean (Phaseolus vulgaris L.). However, its adaptability to different kinds of crops [25] has allowed us to obtain satisfactory results for both calibration and validation datasets, with agreement metrics in line with those reported for other legume species [33,34,35].

Although the extension of the study area was less than 16 thousands square kilometers, the heterogeneity in environmental conditions (soil, climate, and management) led to the identification of more than twenty homogeneous contexts. The variability in the results obtained for both ideotype definition and cultivar recommendations highlighted the marked impact of G × E × M on crop productivity and yield stability. This underlines once more the importance of targeting context-specific ideotypes [15,17] and of accurate characterizations of the environmental context before cultivar choice [26].

Given the complex and dynamic nature of G × E × M interactions—explored here via context-specific ideotyping—the impact of functional traits on the target goal (e.g., maximizing yield) may act counterintuitively and may be specific for the environmental and management conditions considered [55]. In this regard, the comprehensive exploration of the variability in growing conditions in the target area—and of its impact on recommended varieties—allowed us to show that Etna is likely the best choice regardless of soil and climate conditions in case of early sowings, whereas for late sowings there are new varieties (e.g., 18B988 and 17B916) and other cultivars (e.g., Taylor’s and Meccano) that might be a better option for achieving high and stable yields. For instance, in the case of the driest climate (SAM class A) with coarse soil (Figure 6), the slightly earlier flowering and the higher radiation use efficiency of Etna made this cultivar a better choice than Taylor’s, but the situation was overturned in the case of medium and fine soils. This is probably due to the fact that, even if the crop is irrigated, the mild water stress experienced before irrigation is triggered occurs more frequently in coarse soils than in those with higher water holding capacities, which under the dry and warm conditions that characterize the SAM class A (Figure S3) shortened the crop cycle, in turn penalizing crop productivity.

Another example of the marked G × E interactions observed for late sowings is the different ranking achieved by the cultivar Magico across different agro-climatic conditions. This cultivar was indeed among the top-ranked under SAM class A (driest and warmest), but it was clearly penalized under SAM class D. This can be due to the fact that Magico has a long pod filling duration, and the low temperatures that characterize SAM class D contexts (Figure S3) delayed the phenological development and led to maturity occurring under conditions warmer than the crop needs. The negative impact of high temperatures on pod filling due to reduced translocation of photosynthates to reproductive organs [56] is indeed accounted for by the model used in this study [25].

The reliability of the rankings we obtained was supported by the positive feedback received by the stakeholders involved in the study. Cultivar Etna was rated as one of the best in multi-year evaluation trials conducted in the study area, to the point of being used as a benchmark for the evaluation of new lines. As a result, it is currently one of the most widely grown varieties in the study area and, together with Taylor’s and Meccano, it is among the varieties recommended by the Emilia-Romagna Region Agronomic Service. However, preliminary results from field trials carried out in 2017 and 2018 pointed out new promising genotypes, including those identified in this study as being of potential interest (e.g., 18B988).

The methodology developed in this study also indicates second choice varieties (i.e., those with an index of similarity to the context-specific ideotype differing less than 20% from that of the top-ranked variety). This allows supporting—when possible—the recommendation of more than one variety, given that some of the traits that may contribute to drive the farmer’s choice were not included in the analysis. For instance, in case of the borlotto beans, the seed coat colour and other quality traits (e.g., seed shape) are important features for products destined for the dry seed market. Moreover, a larger pool of recommended cultivars would meet the need of sowing more varieties to limit the risk of yield losses due to the susceptibility of a particular genotype to specific biotic or abiotic stressors [26,57]. However, it should be noted that all these traits (qualitative traits and resistance/tolerance traits) can potentially be included in the workflow proposed, the only limit being—as in this study—the availability of functional trait measurements for the varieties of interest.

In particular, the methodology proposed can be easily extended to traits involved with resistance/tolerance to biotic and abiotic stressors. Indeed, examples of ideotyping studies targeting resistance/tolerance traits are already available (e.g., [14,58,59]), and constant efforts are being made by the modelling community to enhance the capability of crop models to reproduce the impact of environmental stressors on crop growth and development (e.g., [60,61]).

5. Conclusions

Variety recommendation is the last step of a breeder’s work and the key starting point of the agricultural production chain. Taking advantage of the capability of crop models to explicitly reproduce G × E × M interactions on large areas and of the available knowledge on the integration of the ideotype concept in cultivar selection frameworks, this study provides the first practical example of exploiting model-based ideotyping for cultivar recommendation purposes.

Results highlighted that even small variation in the environmental and management conditions may have a clear impact on variety suitability. This, in turn, underlines the key role of an accurate characterization of the target pool of environments (climate, soil, management factors) and the importance of grounding the variety recommendation process on context-specific ideotypes.

The crop model STICS was used as a simulation engine, which proved to be able to reliably reproduce bean growth and development despite it being the first time it was used for such a crop. Nevertheless, different crop models can be integrated into the proposed methodological framework. The approach we propose is indeed generic, being easily transferable to other crops, traits, and objective functions (e.g., resource use efficiency, ecosystem services). Given its flexibility and the limited number of inputs required, this approach has the potential for being operationally used to support breeders, seed companies and farmers in the identification of the best varieties for production contexts deriving from high-resolution zonations based on differences in environmental and management factors.