Dissection of genotype × environment interactions for mucilage and seed yield in Plantago species: Application of AMMI and GGE biplot analyses

Genotype × environment interaction (GEI) is an important aspect of both plant breeding and the successful introduction of new cultivars. In the present study, additive main effects and multiplicative interactions (AMMI) and genotype (G) main effects and genotype (G) × environment (E) interaction (GGE) biplot analyses were used to identify stable genotypes and to dissect GEI in Plantago. In total, 10 managed field trials were considered as environments to analyze GEI in thirty genotypes belonging to eight Plantago species. Genotypes were evaluated in a drought stress treatment and in normal irrigation conditions at two locations in Shiraz (Bajgah) for three years (2013-2014- 2015) and Kooshkak (Marvdasht, Fars, Iran) for two years (2014–2015). Three traits, seed yield and mucilage yield and content, were measured at each experimental site and in natural Plantago habitats. AMMI2 biplot analyses identified genotypes from several species with higher stability for seed yield and other genotypes with stable mucilage content and yield. P. lanceolata (G26), P. officinalis (G10), P. ovata (G14), P. ampleexcaulis (G11) and P. major (G4) had higher stability for seed yield. For mucilage yield, G21, G18 and G20 (P. psyllium), G1, G2 and G4 (P. major), G9 and G10 (P. officinalis) and P. lanceolata were identified as stable. G13 (P. ovata), G5 and G6 (P. major) and G30 (P. lagopus) had higher stability for mucilage content. No one genotype was found to have high levels of stability for more than one trait but some species had more than one genotype exhibiting stable trait performance. Based on trait variation, GGE biplot analysis identified two representative environments, one for seed yield and one for mucilage yield and content, with good discriminating ability. The identification of stable genotypes and representative environments should assist the breeding of new Plantago cultivars.


Introduction
Medicinal plants with utilization of a very small cultivation area comprise a huge number of plant species [1,2]. Although only 10% of medicinal species are used commercially, the use of herbal medicines with a diverse array of biological specificities, characteristics and demands is a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 growing [3,4]. This is due to high costs of developing pharmacological products, increasing human interest in natural and organic products, reported side effects of synthetic materials, perceived environmental pollution by the pharmaceutical industry and the difficulties associated with the discovery of new natural medicines [1,[5][6][7][8][9][10]. Plantago ranks sixth in economical trade of medicinal plants [2,11] and has been used traditionally for the treatment of liver disease and cancer, stomach troubles, and for inflammation of the bladder and urinary tract [12]. Its mucilage has a mild laxative effect and it is used commercially to cure diarrhea and colds. The mucilage is also used as an excipient for specific drug formulations and has useful carbohydrate polymers for physiologically active drugs [12,13]. Iran's geographical location is an ideal environment for the growth and cultivation of Plantago. Despite the immense medicinal and export value of Plantago, the productivity of its species is constrained by biotic and abiotic stresses which result in heavy losses in seed/husk quality and yield [4,14]. Moreover, little effort has been devoted to breeding Plantago and attempts to expand its genetic variability, cultivation and domestication have had limited success [4,[15][16][17][18][19][20][21]. Although consumption of herbal medicines is widespread and increasing, domestication and breeding medicinal plants receives little attention. Global climate change has significant effects on environmental conditions and crop production of all crops. Accordingly, environmental factors affecting plant growth and yield should be understood and managed better for more output [22]. Cultivation under diverse environmental conditions opens up the possibility of using breeding programs to solve problems that are inherent in the production of medicinal plants.
Most plant breeding programs are beset by the challenge of genotype × environment interaction (GEI) [23]. GEI is a key factor in the domestication process and cultivation of plants [1,9,[24][25][26][27][28][29]. Approximately 75% of medicinal plant species are collected from the wild and their economic yield is highly affected by environmental conditions [1,3,30]. Seed production, secondary metabolite production and their stability in different climates are target traits in which considerable success can be expected simply by selecting vigorous and stable genotypes, a process that establishes a population adapted to the relevant growing conditions [1,31,32]. The performance of genotypes is always affected by GEI. Several statistical methods have been proposed for analysis of plant stability with the aim of dissecting GEI and stable trait expression across environments. The simple stability coefficient (P i ) shows performance of a genotype in multi-environment trials (MET) [33]. This coefficient defines a superior cultivar as one with a performance near the maximum in MET and measures the rank and superiority of a cultivar, as proposed by Lin and Binns [33,34]. Two frequently used and effective multivariate models for statistical analysis of stability and yield trials are the additive main effects and multiplicative interaction (AMMI) and GGE (genotype (G) main effects and genotype (G) × environment (E) interaction) analysis. These models are used to assess the adaptability and stability of genotypes and to identify mega environments [32,[35][36][37][38][39]. The AMMI model helps to identify genotypes of high productivity adapted to an agronomic zone, with the aim of evaluating environmental [40]. G and GEI are two important sources of variation for evaluating genotypes in MET [31,32,39]. AMMI and GGE biplot analyses combine principal component analysis (PCA) and a graphical explanation of GEI. The AMMI model integrates analysis of variance (ANOVA) and PCA as an efficient procedure for an analysis of the stability of genotypes in a MET [39]. It also quantifies the contribution of each genotype and environment to the total GEI variation [40]. GEI can be better understood by analysis of GGE biplots which in turn facilitates the identification of representative environments, detects the ability of test environments to discriminate and identifies stable genotypes in MET [32,39]. Most studies on medicinal plants have focused on simple analysis of grain yield or secondary metabolites and relatively less effort has been devoted to the advanced analysis of traits in these plants under MET. In addition, few of these analyses have fully dissected GEI or identified representative and discriminating environments to test genotypes [6,[41][42][43][44]. Such information provides a better understanding of the GEI being used in breeding programs of medicinal plants, i.e. Plantago. The objectives of this study were to (1) dissect GEI for grain yield and mucilage performance of Plantago genotypes using AMMI and GGE biplot analyses under various environmental conditions, (2) detect stable genotypes across typical environments and (3) identify representative test environments for future use in breeding programs.

Plant materials
The thirty Plantago genotypes used in this study are listed in Table 1. Seeds of each genotype were provided by the Iran Forests, Range and Watershed Management Organization (IFRWMO).

Field trials and meteorological data
Plantago genotypes were evaluated for mucilage and grain yield under 10 managed field trials. This trait data was compared with trait data collected from Plantago genotypes in their natural habitats. The data for the origin of each Plantago genotype, its code and the Köppen-Geiger climate (KGC) classification [45] of the experimental sites and natural habitats are presented in Table 1. In the KGC classification system, the letters B, C, D, s, w, a, b, k, and h stand for arid, warm temperate, snow, dry summer, dry winter, hot summer, warm summer, cold arid and hot arid climates, respectively. An environment was defined as the combination of location, year and irrigation regime. Managed field trials were performed in two locations, Shiraz and Marvdasht. In Shiraz, Plantago was grown under two irrigation regimes (drought stress and normal irrigation) in three consecutive years (2013-2014-2015) at the Research Farm of School of Agriculture (52˚46 0 N, 29˚50 0 E). A two-year (2014-2015) experiment under two irrigation regimes was performed at Kooshkak Agricultural Research Station (52˚34 0 N, 30˚7 0 E), Marvdasht, Iran. The experimental layout in the managed trials in both Shiraz and Kooshkak was a randomized complete block design (RCBD), each with three replications. A normal irrigation regime had normal levels of irrigation throughout growing season while the drought stress treatment was defined as 50% field capacity (FC) at the 2-3 leafs stage of growth. A seed density of 4 kg ha -1 was used for sowing in soil at a depth of 4 cm. The 1.2 × 1 m experimental plots consisted of two rows 60 cm apart. Seeds were planted in early April in managed field trials in both Shiraz and Kooshkak. 30 kg N ha -1 and 30 kg P 2 O 5 ha -1 were added to the soil at sowing. Weeds were controlled manually after seedling emergence and during the growth cycle. Seed yield (g m -2 ), mucilage yield (g m -2 ) and mucilage content (as percentage per 100 g seeds) were measured for each genotype in each trial.

Mucilage extraction
Mucilage was extracted based on the method of Sharma and Koul [46]. 10ml of HCl (0.1 N) was heated to boiling in a flask and 1 g of Plantago seeds was added to it. Heating was resumed and the process of seed husk dissolution was monitored. When all seeds had changed color, the flask was removed from the heat and the solution was filtered through clean muslin cloth, while still hot. The seeds were washed twice with 5 mL of hot water and the solution was filtered. This process helps separation of traces of mucilage. The dissolved mucilage was mixed with 60 mL of 95% ethyl alcohol, stirred and allowed to stand for 5 h. Finally, the supernatant was decanted and the precipitate was dried in an oven at 50˚C. This represented the total mucilage content. Mucilage yield (g m -2 ) was also calculated following the method of Sharma and Koul [46].

Statistical analysis
Analysis of AMMI model. The AMMI model for the i th genotype in the j th environment is [35,47], where, Y ijr is the mucilage or yield of genotype i in environment j for replicate r, μ is the grand mean, g i is the deviation of genotype i from the grand mean, e j is the environment main effect as deviation from μ, λ k is the singular value for the interaction principal component (IPC) axis k, α ik and γ jk are the genotype and environment IPC scores (i.e. the left and right singular vectors) for axis k, b r (e j ) is the effect of the block r within the environment j, r is the number of blocks, ρ ij is the residual containing all multiplicative terms not included in the model, n is the number of axes or IPC that were retained in the model, and ε ijr is error under independent and identically distribution assumptions, AMMI Stability Value (ASV) was calculated using the formula developed by Purchase et al. [48]: where, SSIPCA1 is sum of squares of interaction principal component analysis 1 (IPCA1) and SSIPCA2 is sum of squares of IPCA2. Sum of the absolute value of the IPC (SIPC) was calculated by a formula developed by Sneller et al. [49], The biplot graph of the AMMI1 (IPCA1 scores vs. additive main effects from genotypes and environments) and AMMI2 (IPCA1 vs. IPCA2) were constructed. P i is the sum of squares of differences of mean genotype i in each environment from the mean of the best genotype in the corresponded environment [33,34], where, X i: is the yield mean of the i th cultivar in the n environments and M is the mean of the maximum response in the n environments. According to Lin and Binns [33], the first part of the Pi expression quantifies the genetic deviation and the second quantifies GEI. Mean rank of each genotype in all environments was calculated as genotype mean rank. GGE biplot analysis. Singular value decomposition (SVD) of the first two principal components was used to fit the GGE biplot model [50], where, Y ij is the trait mean for genotype i in environment j, μ is the grand mean, βj is the main effect of environment j, μ + β j being the mean yield across all genotypes in environment j, λ 1 and λ 2 are the singular values (SV) for the first and second principal components (PC1 and PC2), respectively, ξ i1 and ξ i2 are eigenvectors of genotype i for PC1 and PC2, respectively, η 1 j and η 2 j are eigenvectors of environment j for PC1 and PC2, respectively, ε ij is the residual associated with genotype i in environment j. In GGE biplot analysis, scores of PC1 were plotted against PC2 [32].
The square root transformed form of the data for mucilage content(%) was used in statistical analyses. The AMMI and GGE biplot analyses were performed by GENSTAT software V. 12.0 (GENSTAT 2009) [51]. The mean for environments and genotypes compared using the Least Significant Differences (LSD) test in Statistical Analysis System (SASV 9.2).

Weather conditions and meteorological data at the trial sites and natural habitats
Meteorological data showed that the Plantago genotypes experienced hot temperatures exceeding 30˚C and low rainfall in all experimental trials (Table 2). Two experimental trials at Shiraz received higher annual rainfalls rather than the other 3 trials. Both experimental sites are classified as Bsk in the Köppen-Geiger climate (KGC) classification [45]. Bsk stands for cold arid regions with dry summers. Annual rainfall varied between 214 mm in 2015 at Kooshkak and 373 mm in 2013 at Shiraz. In Plantago, rainfall in the natural habitat regions varied between 123 mm for Minab as a natural habitat of P. ovata and 522 mm for Behshahr as one of habitats of P. major. The meteorological data shows that Plantago species can grow in a wide range of altitudes from below -14 m to above 2000 m.

AMMI analysis of variance
The results of AMMI analysis showed that seed yield, mucilage yield and mucilage content were significantly affected by genotype, environment and GEI ( Table 3). The AMMI analysis of variance indicated that 53.57% of the total sum of squares (SS) for seed yield was captured by the effect of genotype (G) and 30.25% and 14.72% of the total SS were attributable to the environmental (E) effects and GEI, respectively. For mucilage yield, 86.11% of the total sum of squares is justified by G, 5.17% by E, and 7.22% by GEI. G, E and GEI contributed to 91.54%, 1.63% and 4.96% of the total mucilage content variation, respectively. The IPCA9, IPCA8 and IPCA7 were the last significant IPCAs for seed yield, mucilage yield and mucilage content,  (Table 3). AMMI2 (IPCA1 + IPCA2) explained 82.05% of the GEI sum of squares for seed yield, 75.77% for mucilage yield and 53.07% for mucilage content (Table 3).

Discriminating ability and representativeness of environments
The mean and IPCA1 scores of environments and the top rankings genotypes for each environment are displayed in Table 4. The mean comparison of environments showed that drought stress conditions significantly decreased seed yield compared with a normal irrigation regime and with natural habitat conditions. The mean seed yield (36.23 g m -2 ) and mucilage yield (3.96 g m -2 ) were higher in natural habitats of Plantago than means in each of the trial environments (Table 4). Natural habitats yielded higher mucilage content (10.8%) than E6, E7 and E9. However, mean mucilage yield was not significantly decreased under drought stress conditions (E2, E4, E6, E8 and E10). Similarly, mucilage content was not significantly different in drought versus irrigated conditions, although slight increases in content were observed in some drought environments (i.e. E4, E8 and E10) compared with normal irrigation trials (i.e. E3, E7 and E9), respectively, and slight decreases observed in others -E2 and E4, compared with E1 and E3, respectively. In the present study, IPCA1 scores showed that E1, E7 and E6 were main contributors to the stability of genotypes for mucilage yield, seed yield and mucilage content, respectively. Equal zero IPCA1 score indicates highest contribution to genotypes stability, but small contribution to the GE interaction [52,53]. Displacement along the abscissa of AMMI1 biplot graphs indicates main additive effects, whereas displacement along the ordinate shows interaction effects [53]. Environments with IPCA1 scores nearly or equal to zero have small contribution to the interactions and accordingly have large contribution to the stability of genotypes [52,53]. The AMMI1 biplot graph showed that E7 for seed yield, E1 for mucilage yield and E6 and E7, E2, E8 for mucilage content were the largest contributor to the stability of genotypes (Figs 1, 2 and 3; S1, S2 and S3 Tables). The pattern of dispersion and contribution of environments to genotype stability in the AMMI1 biplot graph for mucilage content was completely different from those for seed and mucilage yield (Figs 1, 2 and 3). The AMMI2 biplot graphs show environment scores for seed yield, mucilage yield and content (Figs 4, 5 and 6). In AMMI2, environments which placed near the origin with low scores for IPCA1 and IPCA2 had small contribution to the GE interaction, but large contribution to the stability of genotypes. AMMI2 showed that all environments were far from the biplot origin for seed yield. It also revealed E2, E7 and E9 as relatively high contributor to the stability of genotypes for mucilage content and E4 and E8 for mucilage yield.   In the GGE biplot analysis, the first two PCs explained 93.32%, 99.7% and 96.99% of the total GEI for seed yield, mucilage yield and mucilage content, respectively (Figs 7, 8 and 9; S4 Table). Environments with low IPCA1 and IPCA2 scores which were placed near the origin in the GGE biplot graph have low discriminating ability for genotypes evaluation and high contribution to Genotype × environment interactions in Plantago the stability of genotypes [32,39]. As repeatable interactions divide the target environments into mega environments which discriminate genotypes, genotypes should separately be evaluated for each mega environment [32,39]. In the "which-won-where" GGE biplot graph (Figs 7, 8 and 9), a polygon is drawn by connecting genotypes that are farthest from the biplot origin and all Genotype × environment interactions in Plantago genotypes are surrounded by the polygon. The "which-won-where" GGE biplot graph divided by equality line into sectors in which different mega environments can be detected [31,32]. In the present study, two mega environments were detected for seed yield. The first mega environment included E1, E3, E5 and E7 whilst the second included E2, E4, E6, E8, E9 and E10 (Fig 7).   In the GGE analysis, genotypes were grouped into a single mega environment for mucilage yield (Fig 8). Distances between the environmental vectors show the dissimilarity of environments in discriminating genotypes [32]. Although environments cannot be divided into meaningful mega-environments as they are scattered in one sector of the which-won-where biplot graph, the distance between environments for mucilage yield was based on drought conditions (Fig 8). E2 was identified as a distinct mega-environment for mucilage content in 2013. In 2013, precipitation and mean temperature was different from 2014 and 2015 in Shiraz which contributed to the separation of E2 from other environments in 2013 (Table 2). E5 and E4 were the second mega-environment, and E1, E3, E6, E7, E8, E9 and E10 were the third for mucilage content (Fig 9).
Environment-focused scaling GGE biplot shows average environment axis (AEA) and average environment coordination (AEC). The concentric circles on the comparison biplot graph in Figs 10, 11 and 12 help to visualize the distance of environments to AEA, AEC and the biplot origin. The ideal test environment is represented by the center of concentric circles. For seed yield, E10 had smallest angle with AEA (representative) and was near to the center of concentric circles (Fig 10). The GGE biplot showed that E9 and E10 were the most discriminating and representative environments for the evaluation of genotypes based on mucilage yield ( Fig  11). E8 and E6 were the best representative and discriminating test environments for evaluation of genotypes for mucilage content (Fig 12).

Stability, adaptability and performance of genotypes
Genotypic means and rank, SIPC and cultivar superiority (Pi) are displayed in Table 5. Low SICP reflects the stability of genotypes and low GEI [39]. G26 (P. lanceolata), G16 (P. psyllium) and G1, G2, G4 and G7 (P. major) with low SICP scores were the most stable genotypes for seed yield and G7, G1, G3, G4 and G6 (P. major) for mucilage yield. G4, G5 (P. major), G8 (P. officinalis), G13 (P. ovata) and G18 (P. psyllium) were the most stable genotypes for mucilage content. P i is the mean squares of distance between genotypes i and 'i' where 'i' is the genotype with maximum response over all trials [33]. Genotypes with smaller Pi are closer to the genotype with the maximum yield. Pi represents genotypic superiority in the sense of general adaptability or wide adaptation. The mean rank of each genotype based on P i scores in all environments is presented in Table 5. Based on Pi and mean rank, P. ovata and P. psylilum were the best for seed and mucilage yield, whilst P. ovata and P. lanceolata were the best for mucilage content.

Seed yield (g m -2 )
Mucilage yield (g m -2 ) Mucilage content (%) Ã Differences between genotypes means for each trait (column) that are equal to or less than the LSD 1% is not significant.
Low SICP and ASV reflect the stability of genotypes and low GE interaction. The smaller Pi and mean rank the smaller distance to the genotype with maximum yield and shows best performance. https://doi.org/10.1371/journal.pone.0196095.t005
The "which-won-where" GGE biplot for seed yield (Fig 7) and mucilage yield (Fig 8) placed G12, G13, G14, G15 and G16 (P. ovata) and G17, G18, G19, G26 (P. lanceolata) and G20 (P. psyllium) in the sector of drought environments. These results indicated that such genotypes were well adapted to drought conditions for seed yield. G26 (P. lanceolata) was located near the place of normal environments in the "which-won-where" GGE biplot for mucilage yield (Fig 8). The place of genotypes on the vertices of the polygon in the "which-won-where" biplot shows the best or the poorest performance of each genotype [32]. The perpendicular of the polygon facilitates visual comparison of the distance between genotypes and environments [32]. The place of genotypes and environments in the "which-won-where" GGE biplot for mucilage content differed from those of mucilage and seed yield showing a different response of mucilage content to genotypic and environmental variations. This result agreed with the results of AMMI1 (Fig 3).
The GGE biplot of genotype-focused scaling showed that P. ovata and P. psyllium were at the top of the rankings for seed yield, mucilage yield and content and were closer to the place of the ideal genotype in the biplot graph (Figs 13, 14 and 15; S4 Table). In the graphical analysis of IPCAs, the first principal component (IPCA1) represents cultivar productivity, and the second is associated with cultivar stability [32]. Hence, the GGE biplot showed that the ideal genotype must have a high IPCA1 value (high productivity) and an IPCA2 value close to zero (more stable). The ideal genotype is located in the center of concentric circles of Figs 13, 14 and 15. The genotypes G13 (P. ovata) and G26 (P. lanceolata), for seed yield, G9, G10, G11, G8 (P. officinalis) and G14 (P. ovata) for mucilage yield, and G1, G3, G4, G5 and G6 (P. major), G28 (P. coronopus) for mucilage content were the most stable. The position of these genotypes was in close proximity to the stability axis (average environment axis (AEA)) (Figs 13, 14 and  15). Stable genotypes are desirable only when they have high mean performances and are located closer to the place of ideal genotype in GGE biplot [32].
The response of G5 and G7 (P. major), and G26 (P. lanceolata) to GEI for seed yield (Figs 7, 10 and 13) and G26 (P. lanceolata) for seed and mucilage yield (Figs 8, 11 and 14) reflects the existence of the interspecific variation being used in Plantago breeding.

Discussion
In the present study, the response of thirty Plantago genotypes to environmental conditions was investigated by the AMMI and GGE biplot models and the simple statistic of Pi on the basis of variations in mucilage content and seed yield traits. In the present study, significant differences were found between Plantago species for genotypic means and ranks, SIPC and cultivar superiority (Pi) for mucilage and seed yield. Variation in the interspecific P i index shows the existence of genetic diversity that can be exploited in breeding programs of Plantago. Genetic variation is the fundamental basis for the improvement of Plantago [1]. Results by Canter et al. [1] indicated exploitation of the genetic potential of medicinal plants is still in its initial stages but variation exists between medicinal plants for breeding. The results of the present study revealed that mucilage content was not significantly different in drought versus irrigated conditions, although slight increases in content were observed in some drought environments. It has been reported that drought enhances secondary metabolites due to a stressrelated decline in biomass production associated with an unchanged biosynthesis rate of natural products or an authentic enhancement of the total secondary metabolite content [10,15,[17][18][19][20]54]. The results showed that the effect of genotype on mucilage yield and content was  [45] letter: B: arid, s: summer dry, k: cold arid. more salient than on seed yield whilst seed yield slightly had more response to environmental variations.
AMMI analysis of variance indicated that genotype was the major contributor to the total variation in mucilage and seed yield in comparison to environment and genotype × environment interactions. In a study on rice, genotype accounted for 67.11% of the total variation of grain yield [53]. Similarly, Islam et al. [55] reported that 57.34% of the total sum of squares rice grain yield was attributed to the effect of genotype. Contrary to these results, analysis of yield stability in yellow passion fruit varieties indicated that 61% of the total sum of squares was explained by the environment and 5% and 34% were attributable to the genotype and GEI effects, respectively [52]. One may propose a "70-20-10" saying that a median yield trial has about 70% of variation in E, 20% in GEI and 10% in G, but as the genotypes (in contrast to similar cultivars) and environments become more diverse the share of G and GEI tends to increase [56,57]. Accordingly, for the purpose of plant selection, the rankings of genotypes matter, which are determined by G and GEI and in such cases, reducing or eliminating the influence of environmental main effects is interested [35,58]. Our study indicated that the contribution of AMMI2 to GEI sum of squares was in agreement with the share of AMMI2 in total variance in Jamshidmoghaddam and Pourdad [39] and with the results of a study by Oliveira et al. [55] who showed that AMMI2 biplot may be more accurate to extract GEI variation given as it contains information of two IPCAs and greater pattern proportion compared to the AMMI1 which considers only the IPCA1. The AMMI2 model was simple, allowing conclusions to be made about stability, genotypic performance, genetic divergence between cultivars, and the environments that optimize cultivar performance [59]. In the AMMI2 biplot graph, close genotypes and environments have positive associations and the place of stable genotypes is near the origin of biplot [53,59,60]. In the present study, AMMI2 showed some of genotypes (i.e. G10, G11 and G26) were adapted to drought stress conditions in view of mucilage yield or content whilst they had variations for seed yield. This result shows different response of mucilage as a secondary metabolite versus seed yield to drought stress conditions [61].
The perpendicular of the polygon in the "which-won-where" GGE biplot facilitates the visual comparison of distance between genotypes and environments and helps identify the representativeness of environments and their discriminating ability [32,33]. In the present study, seed yield was mostly affected by the irrigation regime and studied environments were divided into two meaningful mega-environments on the basis of drought conditions. In a study on eighteen wheat genotypes, nine environments were classified into four clusters based on discriminating ability of genotypes and representativeness of environments [33]. Akter et al. [53] identified the most discriminating environments for grain yield using the AMMI analysis of twelve rice genotypes and five environments in Bangladesh. Identifying mega-environments helps evaluate the discriminating ability and representativeness of the environments with a view to detect locations that can be culled without losing important information about the genotypes [62]. In this study, E9 (normal irrigation condition in Kooshkak) was grouped with drought condition trials. This was possibly because of the low precipitation in 2014 in Kooshkak ( Table 2). The information on mega-environments allows breeders to identify discriminating and representative environments that are good test environments for detection of generally adapted genotypes or breeding for adaptation to specific environmental factors [22,62]. In addition, by adding mega-environment boundaries, breeders can determine whether a test location is predictive for a given environment or else is frequently crossing mega-environment boundaries from year to year [63]. The GGE biplot analysis for mucilage content showed that this trait was discriminated approximately by the year effect. This result showed that environmental conditions cause noticeable effects on the life cycles of Plantago that was in agreement with the results of previous studies on medicinal plants [16,17,25,64].
In the GGE biplot graphs, the smallest angle with AEA represents the most representative environment [32]. The most informative and discriminating environment is the farthest from the origin of biplot, whereas a non-discriminating environment which plotted near the origin is not informative, provides little information and should not be used as a test environment in plant breeding programs aimed at cultivar release [32,39]. Discriminating and non-representative environments are candidate test environments to identify genotypes with special adaptability, whilst discriminating and representative environments are good test environments to verify broadly adapted genotypes. Environment-focused scaling GGE biplot revealed that E10 for seed yield and E9 and E10 for mucilage yield were the most discriminating and representative environments. A test environment can be characterized by its similarity with other environments and its discrimination power [65]. The results in the present study showed that drought condition was discriminating for specially adapted genotypes based on seed yield although it was not representative. The results of GGE biplot analysis agreed with previous reports on mucilage content showing it was more responsive to environmental conditions in comparison to seed yield in Plantago [1,25,53,66].
The AMMI model has been shown to be effective as it contributed a large portion of the GEI sum of squares and separating the main and interaction effects. The GGE biplot model is effective for identifying the stable cultivars across environments and identifying the best cultivars for mega-environment differentiation [63,67]. The results showed that the AMMI and GGE biplot models had similar results in view of specific adaptability to environmental conditions. Nevertheless, contrary results were obtained for environmental contribution to the stability of genotypes. For instance, AMMI1 introduced E7 with low contribution to the GE interaction for seed yield and E1 for mucilage yield, but in the AMMI2 and GGE biplot analyses E7 and E1 were high interactive for seed and mucilage yield, respectively. In AMMI2, all environments had high contribution to the GE interaction for seed yield, whilst in GGE biplot E3, E1, E7 and E5 were found as relatively low discriminating with high contribution to genotype stability for seed yield. In AMMI1, stable genotypes were detected only based on IPCA1 scores but relatively different results were found in the AMMI2 and GGE biplot models due to contribution of two IPCAs information in detection of stable genotypes [52,53,55]. On the other hand, the different statistical basis of IPCA in AMMI2 and PC in GGE biplot leads to some differences in stability results [59]. Although differences in the biplot view of genotypic stability were found, the results of AMMI1, AMMI2 and GGE biplot analyses agreed for some of genotypes. For instance, G26 (P. lanceolata) was detected as stable genotype for seed yield based on AMMI1, AMMI2 and GGE biplot. P. ovata, P. psyllium, P. lagopus and P. coronopus were identified as drought-adapted species in AMMI and GGE biplot analyses. The AMMI1, AMMI2 and GGE biplot models also showed similar results for adaptability of genotypes based on mucilage yield and content. Furthermore, the position of G5 (P. major) was relatively similar in AMMI1, AMMI2 and GGE biplot graphs. Also, the distribution of genotypes in both AMMI graphs and GGE biplots with respect to the traits tested demonstrated that genotypes originated from habitats having similar environmental conditions placed close each other. This may demonstrate genotypic characteristics typical for their habitats. This might be due to the effects of habitation variations of a trait in a species which was in accord with the results of Wolff and Delden [68] who reported variations in Plantago lanceolata in the natural habitats was also observed in the greenhouse and experimental garden. Results of the same study revealed that inhomogeneous habitat characteristics and level of phenotypic plasticity affects genetic variation and species microevolution. In a study on wheat, with respect to grain yield, the genotype rankings in the GGE biplot was significantly correlated with the rankings identified in AMMI analysis [46]. In the present study, the results of P i statistic agreed with the results of the AMMI1, AMMI2 and GGE biplot models about top ranked-genotypes in all environments.

Conclusion
Mean seed yield was highly affected by genotype variation and environmental conditions, while mucilage as a secondary metabolite was mostly affected by the effect of genotype. AMMI stability analysis showed that environments divided into normal and drought conditions for mean seed yield. In the AMMI1 biplot, distribution of environments for mucilage content was mostly based on location and the year of experiment. Environments were divided into two meaningful mega-environments for seed yield based on drought conditions. Drought cannot discriminate genotypes for mucilage yield and content, therefore other environmental conditions may be responsible for this purpose. E10 for seed yield and E9 for mucilage yield and content were the best representative environments with highest ability discriminating of genotypes. P. ovata and P. pysllium with the highest seed yield showed positive interaction with drought conditions for seed yield and therefore are good candidates for cultivation and domestication under arid climate of Iran. P. coronopus and P. lagopus also showed positive interaction with drought stress conditions, but they had moderately low mean seed yield. P. major, P. officinalis, P. amplexicaulis and P. lanceolata showed positive interaction with normal irrigation conditions and negative interaction with drought conditions for mean seed yield. The origin of the genotypes G5 and G7 (P. major), G14 and G15 (P. ovata), the identified stable genotypes, belongs to the Bsk climate (arid, summer dry and cold arid) which is similar to the climate classification of the experimental field in the present study. G26 (P. lanceolata) was another useful genotype with respect to the traits tested in this study. This genotype was collected from a climate (the Bsh: arid, summer dry and hot arid) being relatively similar to the climate of the experimental field. The meteorological data indicated that Plantago could grow in wide geographical ranges differing in altitudes, annual rainfall and temperature. These show the importance of climate in cultivation and domestication of Plantago. Results showed that adaptability of some genotypes for mucilage yield and content was different from seed yield. This result emphasizes different response of secondary metabolite and seed yield to stress conditions. Overall, sufficient inter and intra specific variations were found between Plantago species laying the foundation for plant selection and improvement of seed yield, mucilage yield and content in breeding programs aimed at production of new cultivars.