Adaptation of Grass Pea (Lathyrus sativus) to Mediterranean Environments

Grass pea (Lathyrus sativus) is an annual legume crop widely cultivated in South Asia and Sub-Saharan Africa, but in regression in Mediterranean region. Its rusticity and nutritious value is calling back attention for its reintroduction into Mediterranean rain-fed farming systems. We studied the adaptation of a range of breeding lines in multi-environment field testing in Spain and Egypt, showing wide variation for grain yield. Broomrape (Orobanche crenata) infection appeared as the major limiting factor in both countries. Level of broomrape infection was highly influenced by environmental conditions, being favored by moderate temperatures at crop flowering and rain and humidity after flowering. The additive main effects and multiplicative interaction (AMMI) analysis was applied to understand the interaction between genotype (G) and environment (E) on grain yield and on broomrape infection. AMMI analyses revealed significant G and E effects as well as G*E interaction with respect to both traits. The AMMI analysis of variance (ANOVA) revealed that both, yield and broomrape infection were dominated by the environment main effect. AMMI1 biplot for grain yield revealed Ls10 and Ls11 as the accession with highest yields, closely followed by Ls16, Ls18 and Ls19. However, these accessions showed also lower stability, being particularly adapted to Delta Nile conditions. On the contrary, accessions Ls12 and Ls14 were more adapted to rain fed Spanish conditions. Accessions Ls7, Ls1 and Ls3 were the most stable over environments for grain yield.


Introduction
Grass pea (Lathyrus sativus L.) is an annual cool season legume grown for human food but also for animal feed and as forage since ancient times in many of the harshest agro-environments of the world [1]. It is still grown in 1.5 Mha, mainly in South Asia and Sub-Saharan Africa, but is in regression in Mediterranean environments. As world demand for legume protein increases, the rusticity of grass pea recalls its potential for cropping systems in marginal environments as those in Mediterranean cropping systems [2][3][4][5][6]. This called for a renewed interest to gather and characterize the adaptation of landraces and to submit them to breeding in order to exploit the potential of the species [7][8][9][10][11].
However, grass pea suffers from a reputation of being toxic, because under certain circumstances its overconsumption has caused neurolathyrism, a neurodegenerative disease in humans and domestic animals, due to its content of the neuroexcitatory β-N-oxalyl-l-α,β-diaminopropionic acid (β-ODAP) [6,12]. This directed most efforts in breeding to the development of cultivars with a low β-ODAP [1,13]. However, the long-term results of these efforts are often questioned [3,4,6,12]. The β-ODAP content is highly influenced by climatic and edaphic conditions and has a strong Table 1. Description of the environments (combination of location and season) of the trials of the multi-environment study. Summary climatic data corresponding to each growing season are provided (more detailed data used in Figures 3 and 4 are not listed here).

Statistical Analysis
Data for each trait (grain yield and number of broomrapes per plant) were submitted to a combined analysis of variance (ANOVA) with genotype and location-year (environment) as fixed factors using SAS ® 9.3 (SAS Institute Inc., Cary, NC, USA). F-ratios were used to test effects for randomized complete block experimental design [20]. Prior to each ANOVA, tests for normality and equality of variance were conducted for each dependent variable. ANOVA analysis to quantify the contribution of genotype, environment and their interaction (G*E) was performed according to the statistical model: where y egk was the yield or the number of broomrapes per plant of the k block of the g genotype in the e environment; µ the overall mean of the group of genotypes; E e the effect of the e environment, G g the effect contributed by the g genotype; G*E eg the effect interaction between the g genotype and the e environment; B k(e) is the replication effect for k at environment e and ε egk was the residual error effect.

AMMI Analysis
AMMI statistical model structure and computational methods were used in the current study as previously described [21][22][23]. The AMMI method combines the traditional analysis of variance for the main effects of G and E as well as the principal component analysis (PCA) with multiplicative parameters which was used for the partition of the G*E interaction into several components [21]. Analyses were made with the SAS ® 9.4 (SAS Institute Inc., Cary, NC, USA) [24] to graph AMMI biplots.
The AMMI model [21] is given by: where Y ge is either the yield or the number of broomrapes per plant of the g-th genotype in the e-th environment; µ is the grand mean; α g and β e are the genotype and environment deviations from the grand mean, respectively; λ n is the eigenvalue of the PCA axis n; γ gn and δ en are the genotype and environment principal component scores for axis n; N is the number of principal components retained in the model; and ε ge is the residual term. A simple protocol for using AMMI effectively was used here according to Gauch [22]. This protocol includes four steps: (1) analysis of variance, (2) model diagnosis, (3) mega-environment delineation, and (4) agricultural recommendations.
(1) Analysis of variance: the sums of squares (SS) for genotypes (G), G*E signal (GE s ), and G*E noise (GE N ) provide a preliminary indication whether AMMI analysis will be worthwhile. AMMI analysis is appropriate for datasets having substantial G and substantial GE s [22]. To estimate the SS for GE N we multiplied the error mean square from replication by the number of degrees of freedom (df) for G*E. Also we obtained GE S by subtracting GE N from G*E. According to the high SS values of G and GE S , the AMMI model was employed to study the GEI of grass pea yield and number of broomrapes per plant.
(2) Model diagnosis: AMMI constitutes a model family, not a single model. Consequently, model diagnosis was required to determine which member of this model family was the best for a given dataset and research purpose [22]. Malik et al. [25] proposed resampling-based methods for multi-environment trial data with replicates, which are free from distributional assumptions. The R code for the proposed method [25] was used. The p values of the resampling-based tests were derived using B = 5000 bootstrap samples.
(3) Mega-environment delineation: we obtained a ranking table showing the ranks for the accessions in each environment. Mega-environments can be displayed by this table [22]. In order to characterize testing sites in terms of environmental factors which would possibly support the AMMI analysis results, a Non-Metric Multidimensional Scaling (NMDS) ordination [26] was performed. Climate data for each location [maximum, minimum and average temperature, maximum, minimum and average humidity, evapotranspiration and accumulated rain during pre-flowering, at flowering and post-flowering period, also during the growing season] were obtained from the Red Andaluza de Estaciones Meteorológicas (https://www.juntadeandalucia.es/agriculturaypesca/ ifapa/riaweb/web/estaciones) for Spanish locations, and from Metrological Station, Rice Research and Training Center, ARC, Sakha, for Egyptian ones.
(4) Drawing recommendations on best accessions: the goal of yield-trial research is the identification of the best performing genotypes [22]. For this purpose, we used AMMI1 biplot graphs, in which the displacements along the abscissa indicate differences in main (additive) effects whereas; displacements along the ordinate indicate differences in the interaction effects. This was complemented with the rankings of the best accessions per environment in terms of yield and broomrape infection estimated by AMMI1.

Results
Results showed the impact of the environment on the yield and broomrape infection on grass pea accessions. Grain yields of the accessions for each environment are shown in Table 2. Global average for grain yield over environments and accessions was 1883 kg ha −1 . Average yield over environments was highest for Ls10 (2762 kg ha −1 ) and lowest for Ls12 (1368 kg ha −1 ). Average yield over accessions was highest at TOM09 (4103 kg ha −1 ) and lowest at CAMP09 (407 kg ha −1 ). No significant levels of infection of insect pests and fungal diseases were noticed at any site. However, high infection by broomrape infection was recorded in most sites (Table 3). Global average of broomrape infection over accessions and environments was 0.98 broomrapes per grass pea plant. Average infection over accessions was lowest at TOM09 and highest at KAFR08 (0.04 and 1.91 broomrapes per plant, respectively). Average infection over environments was lowest on Ls3 and highest on Ls6 (0.31 and 1.69 broomrapes per plant, respectively). The combined analysis of variance for grain yield and number of broomrapes per plantrevealed that all main effects (environments (E), genotypes (G) and G*E interaction) were statistically significant (Tables 4 and 5). Environment was the most important factor for both traits (Tables 4 and 5), retaining more than a half of the variance contribution while, G and G*E showed similar contribution. Environment explained for 76% of the total variation (G + E + G*E sum of squares (SS)) for yield, whereas G and G*E accounted for 9% and 15%, respectively (Table 4). For number of broomrapes per plant, E, G and G*E accounted for 60%, 17% and 23%, respectively (Table 5).  Tables 4 and 5 also show the degrees of freedom (DF) for the two first principal components of interaction (IPCAs) assigned according to the method of Gollob [27]. In Table 4, the SS for G is 68,746,582, for GE N = 108 × 476,515 = 51,463,620, and for GE S = 113,679,415 − 51,463,620 = 62,215,795. In Table 5, the SS for G is 48, for GE N = 108 × 0.15 = 16.2, and for GE S = 64 − 16.2 = 47.8. In both cases, the GE S was as large as the SS of G. Therefore, AMMI analysis was quite likely to be used worthwhile [13].

Model Diagnosis
According to Malik et al. [25] it is crucial to decide how many significant multiplicative interaction terms to retain. They propose resampling-based methods for multi-environment trial data with replicates, which are free from the distributional assumptions. The p values of the resampling-based tests were derived using B = 5000 bootstrap samples. For yield and K + 1 = 1 the test statistic was of 0.2296693 and a p-value of 0.0000 whereas, for k + 1 = 2 the test statistic was of 0.04037152 and a p-value of 0.5420. For the number of broomrapes per plant and K + 1 = 1 the test statistic was of 0.4515335 and a p-value of 0.0000 whereas, for k + 1 = 2 the test statistic was of 0.04538442 and a p-value of 0.490. For both traits, the resampling-based methods resulted in the same number of significant multiplicative terms. This test suggested that an AMMI model with only one multiplicative term was appropriate for the grass pea dataset. Therefore, the most precise model was the AMMI1.
Correlation coefficients can be calculated between IPCA scores and various genetic or environmental measurements to discern causal factors. For grain yield environmental IPCA1 was positive significantly correlated with evapotranspiration and minimum humidity (r 2 = 0.87, p < 0.01 and r 2 = 0.93, p < 0.01, respectively), and maximum temperature in pre-flowering (r 2 = 0.90, p < 0.01). It was also negative significantly correlated to the maximum humidity in growing season (r 2 = −0.83, p < 0.01). For the number of broomrapes per plant environmental IPCA1 was positive significantly correlated to the growing season evapotranspiration and minimum humidity (r 2 = 0.97, p < 0.01 and r 2 = 0.88, p < 0.01, respectively), and maximum temperature in pre-flowering (r 2 = 0.89, p < 0.01) and it also was negative significantly correlated to the maximum humidity in growing season (r 2 = −0.88, p < 0.01) and minimum temperature in flowering (r 2 = −0.82, p < 0.01). These results are in agreement with NMDS (Figures 1 and 2) and with the Principal Components analysis with climatic parameters (not shown).

Mega-Environments (MEs) Delineation and Description
In the next step of our analysis, we employed the most precise AMMI model to determine the existence of mega-environments [22]. Table 6 is a ranking table listing the top 6 accessions in each of the 7 environments according to AMMI1. This model is of particular interest because is suitable for mega-environment delineation [22]. In our case, AMMI1 is also the model which optimizes predictive

Mega-Environments (MEs) Delineation and Description
In the next step of our analysis, we employed the most precise AMMI model to determine the existence of mega-environments [22]. Table 6 is a ranking table listing the top 6 accessions in each of the 7 environments according to AMMI1. This model is of particular interest because is suitable for mega-environment delineation [22]. In our case, AMMI1 is also the model which optimizes predictive

Mega-Environments (MEs) Delineation and Description
In the next step of our analysis, we employed the most precise AMMI model to determine the existence of mega-environments [22]. Table 6 is a ranking table listing the top 6 accessions in each of the 7 environments according to AMMI1. This model is of particular interest because is suitable for mega-environment delineation [22]. In our case, AMMI1 is also the model which optimizes predictive accuracy according to the resampling-based method [25]. According to AMMI1 analysis, accessions Ls10 and Ls11 were ranked winner for grain yield in most locations (Group 1, sites suffering from broomrape infestation), whereas accessions Ls8 and Ls13 were the winners for this trait in TOM09 (Group 2, the only site with little broomrape). Different cultivars were winners in the case of the number of broomrapes per plant, being Ls3 and Ls1 in the first rank positions (lowest infection) for most environments, with no groups detected.

Selection or Recommendation of the Best Accessions
Accessions falling almost on a horizontal line in the AMMI1 biplot ( Figure 3) have similar interaction pattern lines or environments. For instance, accessions Ls8 and Ls9 were clearly observed from the biplot for similar interaction patterns. Similarly, environments KAFR08 and KFR09 showed similar interaction patterns. On the other hand, accession or environments appearing almost on the vertical line have similar means. For instance, accessions Ls19, Ls18 and Ls16 are shown by the biplot having similar means above general mean of the trial. Also, accession or environment on the right side of the midpoint on the vertical lines have higher mean yield than those on the left side. As a result, it is quite evident from the biplot (Figure 3) that accessions with highest yield were Ls10 and Ls11, followed close by Ls16, Ls18 and Ls19. In contrast, Ls12 and Ls14 were those with lowest yield. Environments on the right hand side of the midpoint of the main effect axis suggested that these environments were more favorable for yield (TOM09, KAFR08 and KAFR09), whereas those to the left (CAMP08, ESC08 and CORD09) were less favorable. Likewise, lines or environments with IPCA1 scores of zero or near zero (either positive or negative) have small interactions and show greater  Environments on the right hand side of the midpoint of the main effect axis suggested that these environments were more favorable for yield (TOM09, KAFR08 and KAFR09), whereas those to the left (CAMP08, ESC08 and CORD09) were less favorable. Likewise, lines or environments with IPCA1 scores of zero or near zero (either positive or negative) have small interactions and show greater stability and adaptation over the test environment. On the contrary, a large IPCA1 scores indicate high interaction effect and reflect adaptation to specific environments. In this present case, ESC08 was the environment with lowest IPCA1 scores whereas KAFR08, KAFR09 and TOM09 were those with highest scores, therefore there were accessions specifically adapted to these environments.
Accessions Ls7, Ls1 and Ls3 recorded the lowest IPCA1 scores and showed small interactions (close to vertical axis). On the contrary, accessions Ls18, Ls16, Ls11 and Ls10 recorded the highest IPCA1 score. This implied that high interaction effect of lines had specific stability and adaptations to particular environments.
The stability of responses to broomrape can be seen in the AMMI1 biplot ( Figure 4). Accessions Ls6, Ls19, Ls12, Ls17 and Ls14 interacted positively with KAFR08 and KAFR09 (all in the right upper corner of the biplot), indicating that they are best suited for these environments, but negatively with the remaining environments. Similarly, accessions Ls3, Ls1, Ls18, Ls7 and Ls10 interacted positively with CAMP08, ESC08, ESC09, and CORD09 (all in the left down corner), but negatively with KAFR08 and KAFR09. The group of lines Ls11, Ls15, Ls4 and Ls5 and the group of lines Ls9, Ls16, Ls2 and Ls13 have low interaction and hence stable over environments.

Discussion
Common experience in Mediterranean countries is a severe decrease in grass pea cultivation during the last 50 years, with little or negligible trade. However, farmers share an appreciation on derived food products and show interest in the crop, provided that better adapted cultivars are made available. As a highly nutritive food/feed crop and its high drought tolerance with minimal input requirements for its cultivation, grass pea provides food and nutrition security to many low-income communities.
Great morphological variation is reported in grass pea, showing a clear grouping in two major types, one covering blue flower and colored seeds, and other white flowers and seeds [28]. The first  With the help of the rankings shown in Table 6 we can recommend for broomrape infection Ls3, Ls1 and Ls18 in spite of their low value in Figure 4, showing that they are better adapted to Egyptian than to Spanish locations. The most interesting accessions would be those with high yield and low broomrape infection, particularly Ls10, Ls11 and Ls18.

Discussion
Common experience in Mediterranean countries is a severe decrease in grass pea cultivation during the last 50 years, with little or negligible trade. However, farmers share an appreciation on derived food products and show interest in the crop, provided that better adapted cultivars are made available. As a highly nutritive food/feed crop and its high drought tolerance with minimal input requirements for its cultivation, grass pea provides food and nutrition security to many low-income communities.
Great morphological variation is reported in grass pea, showing a clear grouping in two major types, one covering blue flower and colored seeds, and other white flowers and seeds [28]. The first group is typical of Indian subcontinent whereas the second has more western distribution. Our Ls accessions correspond to the first group, whereas Ls19 has larger and whiter seeds, the type preferred by Mediterranean consumers. Mediterranean accessions are regarded also to have lower content of β-ODAP [13], although we did not assess this trait. The β-ODAP content is highly influenced by climatic and edaphic conditions and has a strong G*E effect [3,14,15]. This, together with the current understanding that the disease is caused only when there is overconsumption of grass pea in a non-balanced diet, made breeders reconsider yield stability as a higher concern. We therefore focused in this study on adaptation of current breeding lines to different farming systems in the Mediterranean Basin, identifying the more productive for each of the studies locations and the more stable ones. AMMI1 biplot (Figure 3) allowed identification of the accessions better adapted to Delta Nile conditions (Ls10, Ls11, Ls16, Ls18 and Ls19) and to rain fed Spanish conditions (Ls12 and Ls14) as well as those most stable over environments (Ls7, Ls1 and Ls3).
We also paid attention to other factors potentially affecting yield, such as pests and diseases. Grass pea has been reported to be affected by a number of pests and diseases, including aphid, weevil, rust, powdery mildew or fusarium wilt [6,29]. However, in Mediterranean regions, the weedy root parasite crenate broomrape is regarded as a major constraint for grass pea cultivation [30]. Our current results confirm this, as we did not identify any other pest or disease at worrying levels in any of the locations. Grass pea appears as very sensitive crop to broomrape infection, suffering higher yield penalty than pea or faba bean [31], with 100% reduction at severities higher than four parasites per plant. In addition to the host, infection severity of broomrape strongly depends on parasitic seedbank density and on environmental factors such as temperature and rain. Level of broomrape infection was favored by moderate temperatures and high humidity at crop flowering (Figure 2), what coincides with broomrape seed germination and establishment, and by higher pluviometry after crop flowering, what coincides with broomrape development and emergence. Temperature effects have already been reported, with cool winters acknowledged to reduce infection [32][33][34]. Water availability is indeed acknowledged as an important factor for broomrape development [32,33,35]. Small levels of resistance have been reported in grass pea [20] which would be the result of a combination of different escape and resistance mechanisms, acting each one at different phases of the infection process as is the case in most crops [36]. AMMI1 biplot (Figure 4) allowed the identification of accessions performing better in Egypt (Ls6, Ls12, Ls14, Ls17and Ls19), those performing better in Spain (Ls1, Ls3, Ls7, Ls10 and Ls18) and those stable over environments ((Ls2, Ls4, Ls5, Ls9, Ls11, Ls13, Ls15 and Ls16). The most interesting accessions would be those with high yield and low broomrape infection (Table 6), particularly Ls10, Ls11 and Ls18. Availability of accessions with high and stable yield offers hope for the reintroduction of the crop in Mediterranean conditions. Author Contributions: D.R., A.A.E. and late A. Moral designed and performed the trials; D.R. and F.F. analyzed the data and wrote the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the Spanish Agencia Estatal de Investigación (AEI) grant AGL2017-82019 to finalize the data analysis.