Evaluation of yield, yield stability, and yield–protein relationship in 17 commercial faba bean cultivars

Faba bean is a legume crop with high protein content and considerable potential for wider cultivation in cool climates. However, it has a reputation for having unstable yield with large interannual variability, mostly attributed to yearly variation in rainfall. In this study, 17 commercial cultivars of faba bean were evaluated for seed yield, yield stability and the relationship between seed yield and protein content at four locations in Denmark and Finland during 2016–2018. We found that location and year effects accounted for 89% of the total seed yield variation. Cultivar × environment (GxE) interactions were small (2.4%) and did not cause reranking of cultivars across environments. Yield stability contributed little to the mean yield of the cultivars, as high‐yielding cultivars consistently outperformed the lower yielding genotypes, even under the most adverse conditions. Similarly, GxE effects on protein content were limited, and we found an overall negative correlation of −0.61 between seed yield and protein content for the cultivars and environments studied. These data may be helpful for selecting cultivars for field use or for use in breeding programmes, considering that future faba bean pricing could depend on both protein quantity and concentration.

Faba bean is cold-hardy and grows best under cool and moist conditions, making it a very attractive legume crop for cool climates where soybean does not thrive (Duc et al., 2015). It has shown good potential for use in pig feed (Partanen, Siljander-Rasi, & Alaviuhkola, 2006), and its average protein content is higher than that of other common food legumes such as pea (Feedipedia.org, 2020).
The protein content has been found to range from 22% to 38%, with spring cultivars generally showing higher levels than winter cultivars (Bond, 1977;Griffiths & Lawes, 1978). Replacing part of the soybean import in European countries with locally grown faba bean will benefit not only the environment but also the national economies (Stoddard et al., 2009).
However, the yield of faba bean, like that of many other legumes, is considered unstable because of large interannual variability in yield as compared with nonlegume species (Cernay, Ben-Ari, Pelzer, Meynard, & Makowski, 2015). Much of this instability is attributable to the contrast between spring sowing of most legume crops and autumn sowing of most of the cereals with which they are compared . Faba bean is considered drought susceptible, and variation in the amount of rainfall has been reported as a major cause of seed yield instability (Link et al., 1999). To make faba bean a more attractive crop to agronomists and farmers, the identification of cultivars with high seed yield and high yield stability is desirable. When collecting data from multienvironment yield trials, yield stability can be investigated as the ability of some genotypes to perform consistently over a wide range of environments (Finlay & Wilkinson, 1963). The demands for cultivars better adapted to changing conditions is of great interest even in what seems to be consistent environments, as climate conditions change from year to year and are expected to fluctuate more in the future as a consequence of global warming. There is already an indication that grain legume yield variability has increased in the last 60 years (Reckling, Döring, Bergkvist, Chmielewski, et al., 2018).
Several methods may be used to assess yield stability. Static (type 1) yield stability refers to genotypes that perform consistently in all environments, that is, environmental stability, whereas dynamic (type 2) yield stability refers to genotypes that show no genotype × environment (GxE) interaction (Annicchiarico, 2002). Among the most common measures of static yield stability are the coefficient of variation (CV), whereas the Finlay-Wilkinson regression, which works by regressing performance of each genotype on the environmental means, is widely used to describe the GxE interaction. Two kinds of GxEs are widely recognized: one that causes reranking of cultivars between environments and shows a clear pattern of adaptation and another that shows no reranking but where differences between cultivars are larger in some environments than in others. Covering both the static and dynamic types of stability is expected to benefit the identification of cultivars with stable yield.
The objectives of this study were to evaluate the seed yield of 17 commercial cultivars of spring-type faba bean across different Nordic environments, to determine the stability of yield for each cultivar and to assess the relationship between yield and protein content of seeds.

| Plant material and growing conditions
Seventeen commercial spring-type cultivars of faba bean originating from cool climates (Table 1) were grown in four Nordic locations during 2016-2018. Three of the locations were in Denmark: Dyngby (55.942 N,10.208 E),Nørre Aaby (55.470 N,9.845 E) and Horsens (55.830 N,9.950 E). The fourth location was the University of Helsinki research farm at Viikki, Helsinki,Finland (60.224 N,25.021 E).
The combination of a year and location was considered as an environment. The trials were sown in April to May and harvested in August to September. Precise dates of sowing and harvest of each trial can be found in Table S1. No plots were evaluated for yield in Horsens and Nørre Aaby during 2016, owing to a treatment mistake, and no harvest took place in Viikki during 2017, owing to autumn weather damage, so there was a total of nine environments for yield data. Plots were harvested with plot scale that combines harvesters.
All trials were rainfed. Treatments of trials were site specific and included application of fertilizers and pesticides at all locations.
Fungicides were applied to all trials, except those in Dyngby and Viikki during 2016. Comprehensive information about treatments at each location can be found in Table S2. None of the fields showed any signs of uncontrolled strong biotic stress during the growth period.

| Phenotyping
Seed moisture content (%) and protein content (%) were determined with a near-infrared (NIR) sensor (DA 7440, Perten, Stockholm, Sweden; or Infratec 1,241, FOSS, Hillerød, Denmark) on samples from Dyngby in all years and Horsens in 2017. The calibration model used for determination of nitrogen content was artificial neural networks (ANNs). The NIR protein content was verified on a sample size of n = 220 and showed an R 2 of 0.91. The NIR moisture content was verified on a sample size of n = 221 and showed an R 2 of 0.99.
For samples from Finland, protein content was determined using the Kjeldahl N × 6.25 method, and moisture content was determined with a Dickey-John moisture meter (Dickey-John Inc., Auburn, IL, USA). After seed cleaning, plot yields (g/m 2 ) and protein contents were corrected to 15% moisture content. This gave a total yield data set of 487 observations and a protein data set of 270 observations. Additional traits recorded in some of the trials were thousand grain weight (TGW), days from sowing to maturation of plots (plots were considered mature when 95% of the pods were black), days to flowering (considered when 90% of the plot has flowered) and duration of flowering.

| Yield analyses
For the analysis of yield measured across different years and locations, the following linear mixed model was applied using the 'lme4' package (Bates, Mächler, Bolker, & Walker, 2015) in R (v. 3.5.1) (R Core Team, 2018).
where Y ijkn is a random variable representing the seed yield (g/m 2 ), C i is the effect of the ith cultivar (i = 1, 2, …, 17), A j is the effect of the jth year ( j = 1, 2, 3), L k is the effect of the kth location (k = 1, 2, …, 4), AL jk is the effect resulting from year × location interaction, and α i x jk is used to model GxE interaction effects being the product of the cultivarspecific regression coefficient on the environmental mean yields of the ith cultivar (α i ) and the environmental mean seed yield in the environment corresponding to the kth location and jth year (x jk ) as described by the Finlay-Wilkinson regression (Finlay & Wilkinson, 1963). CA ij and CL ik model the effect resulting from cultivar × year interaction and cultivar × location interaction, respectively, meaning that together they model part of GxE interaction that is not captured by the linear α i x jk term. ε ijk is the residual error term of the ith cultivar in the jth year and kth location and μ is the overall mean.
T A B L E 1 Cultivar origin, specification of zero tannin types (0/1 represent presence/absence of tannin), average thousand grain weight (g), days to maturation (number of days from sowing), days to flowering (number of days from sowing), duration of flowering (number of days) and correlations between cultivar means of traits with seed yield (g/m 2 ) and seed protein content (%) The model assumes normality of data and independence of residuals. Visual inspection of QQ-plots and plots checking for homoscedasticity did not give rise to any concerns regarding normality. From the model, restricted maximum likelihood (REML) estimates of different variance components were extracted. The statistical significance of interaction terms modelled as random effects was tested using the analysis of variance (ANOVA) function in R to compare the loglikelihood of full model 1 fitted using ML with a model where the random effects under evaluation were left out one at a time. p values report if keeping the random effect in the model resulted in a significantly improved model in log-likelihood terms. When testing statistical significance of main effects such as year, location and cultivar effects, interaction terms including the effects were removed from the full model before the main effect was dropped and significance was tested. In addition, estimates of variance were used for calculation of differently defined broad-sense heritability terms, including heritability of single plots, heritability of cultivar means and heritability of cultivar means within a single year and location. The equations used to obtain the heritability estimates were the following: where σ 2 c , σ 2 a , σ 2 l , σ 2 e , σ 2 al , σ 2 ca , σ 2 cl and σ 2 αx are the variances of the cultivar, year, location, residual, year with location interaction, cultivar with year interaction, cultivar with location interaction and the cultivar by Cultivar yield means were compared with a Tukey test in R using the built-in Tukey honestly significant difference (HSD) function.

| Evaluation of yield stability
To examine the yield stability of the different cultivars across different environments, we used the CV as a parameter to evaluate type 1 stability and Finlay-Wilkinson regression to evaluate type 2 stability.

Coefficient of variation
The CV measures the variation of a cultivar across environments. It can be defined as follows: where σ i and μ i are the standard deviation and mean of yield data for the ith cultivar, respectively. Cultivars with the highest stability will give rise to the lowest coefficients of variations.

Finlay-Wilkinson regression coefficient
Yield stability was accessed across different environments by the Finlay-Wilkinson regression coefficient (Finlay & Wilkinson, 1963).
The yield response of a given cultivar can be represented as follows: where R ie is the modelled yield response of the ith cultivar in the eth environment, β i is the intercept value of ith cultivar, x e is the mean yield in the eth environment and α i is the regression coefficient of the ith cultivar and is a measure of genotype by environment interaction.
It applies that α i Á x e + β i 6 is equal to C i + α i x jk 1. The lower the regression coefficient, the greater the resistance of the cultivar to environmental change, that is, the higher its stability. Wald tests were applied to calculate if regression coefficients were statistically significant from 1.

| Correlations
Genotypic correlation coefficients between seed yield and protein content of seeds were estimated as the Pearson correlation between the sets of cultivar BLUPs obtained by using protein or yield as the random variable in Model 1. All Pearson correlations and related p values were calculated using the "correlation" function in the "agricolae" R-package (de Mendiburu, 2010).
The R-script used for statistical analyses can be found at https:// github.com/cks2903/Faba_bean_yield_2019 3 | RESULTS

| Seed yield
The 17 commercial spring-type faba bean cultivars (Table 1)  The seed yield showed a significant positive correlation with TGW (0.81) and days to maturity (0.68). The germplasm collection also showed variation in earliness and duration of flowering, but we did not find seed yield to be significantly correlated with these traits (Table 1).
The environments included years with very different weather conditions ( Table 2). The growing season in 2018 proved to be unusually dry and was characterized by higher average temperatures, more hours with sunshine, fewer days with precipitation, lower average amount of rainfall, lower humidity and many days with high risk of drought (drought index > 9) than in 2016 and 2017. Environmental mean yields showed a significant positive correlation with the days with precipitation, the average amount of precipitation and the length of the growing period. Significant negative correlations were found with the average daily temperature, average hours with sun and number of days with drought index larger than 9. We found no significant correlations between weather data and protein content of seeds.
Estimates of cultivar, year, location, year and location interaction, cultivar by Finlay-Wilkinson regression, cultivar × location interaction, cultivar × year interaction, and residual variances and the proportion explained by the different components are presented in Table 3. The seed yields of the cultivars were significantly affected by all factors included in the model. The largest proportion of the variance in seed yield was explained by environmental factors including year (60.7%), location (20.4%) and year × location interaction (7.5%), which together accounted for 88.6% of the total phenotypic variance in yield. GxE interaction effects accounted for 2.4% of the total variance.
The cultivar was estimated to account for 5.3% of the variance.

| Yield stability
The CV% (a measure of type 1 yield stability) varied from 36.9% (Snowdrop) to 45.7% (Gloria) with a mean of 40.8%. Plotting the CV of each cultivar against its mean yield allowed us to identify 749-13, Fuego, Lynx, Pyramid and Taifun as high yielding with low CV% (Figure 1b). Of these, Lynx and Pyramid showed significantly larger yield means than the lowest yielding-cultivar Kontu (Figure 1a). No significant correlation was observed between CV and seed yield.
We found that 50% of all GxE interaction involved in seed yield response, that is, the sum of the variances explained by cultivar × year interaction, cultivar × location interaction and cultivar by FW regression, could be explained by the Finlay-Wilkinson regression coefficients, a measure of type 2 or dynamic yield stability (Table 3).

| Yield-protein relationship
The cultivar averages of protein content ranged from 28.3% in 749-13 to 31.9% in Gloria. The minimum protein content was observed in Lynx (26.6%) and the maximum in Mistral (34.9%). Climatic information on the nine environments and correlations with environmental mean seed yield (g/m 2 ) and protein content of seeds (%) effects for both protein content and seed yield. No significant phenotypic correlation between seed yield and protein content of seeds was observed on a single plot level, but we observed a correlation of −0.68 (p < 0.01) when the two traits were averaged across cultivars (Table 1). In addition, the protein content of seeds averaged across cultivars was significantly correlated with TGW (−0.60, p < 0.05) ( Table 1). The Finlay-Wilkinson regression coefficients showed no significant effect for seed protein content (Table S3), meaning that the type 2 stability of protein yield (protein content of seeds multiplied by seed yield) should be expected to parallel to the stability of seed yield alone.

| DISCUSSION
The germplasm collection included in this study consisted exclusively of commercial spring-type cultivars bred to perform well in cool Nordic climates. For this reason, the results do not represent all variability F I G U R E 2 Regression of mean seed yield (g/m 2 , calculated at 15% H 2 O) of a cultivar in a given environment against the environmental mean yields. The black lines display the average cultivar × environment (GxE) effect. Values in the lower corner refer to the regression coefficient (α) and the R-value (R) of the fitted line (green). * Significant at p < 0.05; ** significant at p < 0.01 F I G U R E 3 Genetic correlation between the seed protein content and yield (both at 15% H 2 O) of the 17 commercial faba bean cultivars. The shaded area specifies the 95% confidence interval of the regression line within the faba bean crop but are specific for this set of cultivars. The cultivars were chosen to reflect several qualities of commercial importance including both good agronomic and disease-related characteristics.
In agreement with our findings, previous studies have found the main factors determining faba bean seed yield to be environmental (Fikere, Tadele, & Letta, 2008;Temesgen, Keneni, Sefera, & Jarso, 2015). The remarkably high proportion of seed yield variance explained by year alone can be attributed to the large annual climatic variation caused by the unfavourable weather conditions during 2018 compared with the other growing seasons. The lowest average daily precipitation, fewest days with precipitation, highest average daily temperatures, hours with sunshine and number of days with high risk of drought occurred during the 2018 growing season, which were connected with lower seed yields than observed during other test years. The presence of these drought indicators of 2018 is expected to cause abiotic stress, which induces physiological and biochemical changes resulting in a lower seed yield (Ammar et al., 2015;Link et al., 1999). Deficiency of water, especially during flowering, early podding and grain filling can cause large reductions in faba bean seed yield (Katerji, Mastrorilli, Lahmer, Maalouf, & Oweis, 2011;Mwanamwenge, Loss, Siddique, & Cocks, 1999). In addition to having an effect on yield, water deficiency tends to increase the protein content in faba bean seeds (Alghamdi, 2009). We also found the seed yield to be highly correlated with the days of growth, which were lower during 2018 as drought stress is known to trigger maturation (Ammar et al., 2015).
The different geographical locations in Denmark had similar annual weather conditions. However, within years, the location in Finland was characterized by lower relative humidity and a higher average temperature and daily amount of precipitation in 2018 compared with the Danish locations. As relative humidity did not show a significant correlation with the seed yield, and temperature and average amount of precipitation seem to exhibit opposing effects on yield, most of the relatively large effect of location is expected to reflect differences in the field conditions across locations rather than geoclimatic differences. In this context, it is worth noting that even in the extremely dry year 2018, the average yield at the location in Dyngby (DK) was 424 g/m 2 , which was only slightly below the overall average (443 g/m 2 ), indicating that faba bean can produce decent yields even under very dry conditions.
In this study, we identified Fuego, Lynx, Pyramid and Taifun as cultivars that were associated with a positive effect on seed yield, while having a high type 1 stability compared with other high-yielding cultivars. None of the four cultivars were less sensitive to GxE interaction than average. For the analysis of GxE interaction or type 2 yield stability, we found a relatively small but statistically significant amount of GxE interactions that did not cause reranking of cultivars in different environments. High-yielding cultivars consistently outperformed the low-yielding cultivars in all environments, and most cultivars showed good adaptation to all environments as indicated by their Finlay-Wilkinson regression coefficients not being significantly different from 1. However, Gracia, Pyramid and Vertigo had better adaptation to high-yielding environments as indicated by regression coefficients above 1 (Finlay & Wilkinson, 1963). In addition, Kontu, Snowdrop and SSNS-1 did not have the ability to exploit the highyielding environments, having Finlay-Wilkinson regression coefficients below one. Cultivars that perform reasonably well under all environments evaluated in this study show the ability to produce reasonable yield in very diverse climatic conditions, so they would be expected to perform well across all relevant Nordic countries where spring-type faba beans are grown.
The relatively low Finlay-Wilkinson regression coefficients and the lack of cultivars that clearly seemed to be worst yielding in one environment and best yielding in another are likely due to the germplasm sample only including spring-type cultivars bred to perform well in the cool-climate locations included in this study and the climatic similarity of the different geographical locations included. It has previously been shown that growing different faba bean cultivars in climatically contrasting sites produce large crossover GxE interactions and call for specific breeding for each geoclimatic area (Annicchiarico & Iannucci, 2008;Flores et al., 2012). In contrast to these studies, we found the effect of genotypes to be larger than the effects of GxE interaction. The proportion of seed yield variance determined by genotypes (5.3%) is the broad-sense heritability estimated based on the replicated single plot data and not the cultivar means. In reality, breeding is done on a genotype level, so the heritability of cultivar means is of greater interest. We estimated the broad-sense heritability of the cultivar mean seed yield to be 16% across environments and 84% within a fixed location-year, where the latter is relevant for comparing cultivars within an environment to determine how much of the seed yield response is due to genetics. The heritability of the cultivar mean yield was found to be relatively low compared with earlier reports on faba bean (0.62-0.77) (Alan & Geren, 2007;Toker, 2004).
The value of the heritability was affected by both the genetic material chosen and the large diversity of environments. The data set used for yield evaluation in this study contained more observations from the dry year 2018 than from 2016 or 2017, and the broad-sense heritability of faba bean yield has been found to be higher under well-watered conditions than under drought conditions (Link et al., 1999).
It has previously been reported that there is a fundamental conflict between breeding cultivars with both early maturation and high seed yield, as the correlation between growing period and yield is typically found to be significantly positive (Prohens, Nuez, & Carena, 2008). Our findings support this relationship, which is clearly demonstrated by the earliest maturing cultivar Kontu having the lowest yield on average. In addition to being associated with higher yield potential, late maturation in faba bean is associated with lower yield stability due to risk of lodging (Döring, 2015).
Among the 17 cultivars included in the study, four were tannin free. Low-tannin cultivars have in general been associated with higher seed protein levels (Crépon et al., 2010). In agreement with this, we found three out of these four to have a protein content above average with especially Mistral and Gloria showing remarkably high levels.
The remaining tannin-free cultivar Taifun was found to behave differently, showing a protein content below average but a seed yield above average. The absence of tannin in faba bean has been found to be caused by either one of two recessive and complementary genes named zt1 and zt2 (Picard, 1976). Carrying the zt2 gene has been reported to have a stronger effect on protein than having the zt1 gene, but as all four tannin-free cultivars studied were believed to carry the zt1 gene, this does not explain the low protein content of Taifun (Duc, Marget, Esnault, Le Guen, & Bastianelli, 1999).
Although it has earlier been reported that many legume crops such as pea and chickpea show a negative correlation between seed yield and protein content of seeds, others have shown an absence of correlation attributed to their capacity for biological nitrogen fixation (Frimpong et al., 2009;Jha, Arganosa, Tar'an, Diederichsen, & Warkentin, 2012;Stoddard, Marshall, & Ali, 1993). Many studies in faba bean have reported a similar absence of correlation (Picard, 1977;El-Sherbeeny & Robertson, 1992). However, some studies show an indication of a negative correlation, although mostly not significant, and other studies have even reported a significant negative correlation (Barłóg, Grzebisz, & Łukowiak, 2019;Bond, 1977;Lizarazo et al., 2015), which is consistent with our findings of a significant negative genetic correlation as well as a phenotypic correlation calculated on cultivar means.
The conflicting reports on the relationship between yield and protein content of seeds suggest that our finding should be considered specific for the set of cultivars studied. We report a negative relationship between cultivar effects of protein and yield responses that is highly dependent on the five cultivars Gloria, Kontu, Mistral, Snowdrop and SSNS-1, which show low yield and high protein content. For this reason, our findings do not provide evidence that the genetic cause of both traits can be traced to the same genetic region or gene as is the case in linkage disequilibrium or pleiotropy. Hence, the information on the yield-protein relationship provided here will be more relevant to farmers than to breeders. It is possible that the negative relationship is an indirect effect of other traits exhibiting opposing effects on protein content and yield. In this context, the correlation between TGW and both traits was found to reflect the negative relationship of yield and protein content in this specific set of commercial cultivars. In particular, the three lowest yielding cultivars Kontu, Snowdrop and SSNS-1 were all found to have remarkably low average TGW values (<400 g), whereas Taifun was the only high-yielding cultivar associated with a TGW below average. In addition, all cultivars except Mistral with average seed protein content above 30% had values of TGW smaller than average. The negative relationship of seed yield and protein was also reflected in the time to maturation, although the relationship with protein content was not found to be significant.
As a result of the narrow set of geoclimatic environments included in this study and the consequently low GxE interaction terms reported, it is not unlikely that some of the cultivars, especially the ones not showing ability to exploit the high-yielding environments (Kontu and Snowdrop), might show better adaptation to other environments. Such differences in cultivar adaptation could also be part of the explanation of the negative yield-protein correlation observed.

Consistent
with previous findings (El-Sherbeeny & Robertson, 1992), we found seed yield and protein yield to be highly correlated 0.93 (p < 0.001), meaning that, judging from this set of faba bean cultivars, growing high-yielding cultivars will generally also result in higher total protein production.

| CONCLUSION
In this study, we observed large variation in seed yield and yield stability of the 17 commercial faba bean cultivars evaluated. Seed yield was found to be a trait highly influenced by environment; that is, the proportions of variance explained by environmental factors and GxE interactions were 89% and 2.4%, respectively. It was possible to identify high-yielding faba bean cultivars with high static yield stability (Fuego, Lynx, Pyramid and Taifun). In addition, GxE interactions caused no reranking of cultivars in different environments, and highyielding cultivars consistently outperformed the lower yielding ones in all environments. In the geoclimatic locations studied, it could therefore be beneficial to grow high-yielding cultivars with a high GxE interaction of seed yield because there is no trade-off with respect to seed yield under suboptimal conditions.

ACKNOWLEDGEMENT
The work was funded by NORFAB: Protein for the Northern Hemisphere, Innovation Fund Denmark Grant 5158-00004B.