Integration of Various Stability Models to Identify High Yielding and Stable Genotypes of Pigeonpea [Cajanus cajan (L.) Millspaugh]

Pigeonpea [Cajanus cajan (L.) Millspaugh] is a climate resilient multipurpose food legume serving as a lifeline to resource poor farmers of India as its cultivation requires minimum external inputs.  It is an ideal legume crop for rainfed tropics and sub-tropics because of its tremendous ability to survive various kinds of biotic and abiotic stresses (Saxena et al., 2018). The major pigeonpea producing countries includes India, Myanmar, Kenya, Tanzania, Malawi, Uganda and Mozambique (FAOSTAT, 2017). India contributed about 72% of global pigeonpea production (Sharma et al., 2019). In India, pigeonpea is mainly consumed as dal in northern and western India, sambar in Southern India and as fresh vegetable in Karnataka and Gujarat state. In India during 2017-18 pigeonpea was annually grown on 4.43 million hectare area with an annual production of 4.25 million tones and a productivity of 960 kg/ha (Anonymous, 2018). The development of novel high yielding varieties well adapted across environments is the ultimate aim of any pigeonpea breeding programme (Meena et al., 2017). Grain yield is very complex trait as it is governed by polygenes and highly influenced by prevailing environmental conditions and in a more precise way it can be concluded that grain yield is governed by genotype (G), environment (E)  and their interaction (G x E) effects (Vaezi et al., 2019). The presence of G x E in multi-environment trials is expressed either as inconsistent responses of different genotypes (relative to others) due to changes in genotypic rank, or as the absolute difference between genotypes without rank change (Crossa, 2012). The G x E effect tends to decreases the relationship between phenotypic and genotypic values and hence hinders the selection of superior varieties (Allard and Bradshaw, 1964). 
       
Genotype x environment interaction complements the selection process and recommendation of a genotype for a target environment and hence a number of different methods have been proposed to estimate G x E interaction (Gauch and Zobel, 1996; Gauch, 2006). These methods can be broadly classified into two main categories i.e. Parametric and Non-parametric methods (Huehn, 1996). Parametric methods mainly includes ecovalence (W²_i) method (Wricke’s, 1962), regression coefficient (b_i) and deviation from regression (S²d_i) method (Eberhart and Russell, 1966), stability variance (σ_i²) method (Shukla’s, 1972), coefficient of variance (CV_i) method (Francis and Kannenberg’s, 1978). An additive main effects and multiplicative interaction (AMMI) model based stability parameters such as ASV (AMMI Stability Value) and YSI (Yield Stability Index) were also of parametric model (Purchase et al., 2000; Bajpai and Prabhakaran, 2000).  The  non-parametric methods includes S(i) measures  of  Huehn’s (1990) and Nassar and Huehn’s (1987), NP⁽ⁱ⁾ measures of Thennarasu’s (1995) and KR or rank-sum of Kang’s (1988).  These traditional parametric as well as non-parametric statistical methods for estimation of G x E interaction are either ANOVA based or PCA based and hence, have some limitations.  As ANOVA is an additive model and can describe only the main effects effectively, whereas PCA only describe multiplicative components (Zobel et al., 1988).  Another problem with these traditional methods is that they are univariate, but the genotypic response to an environment is multivariate (Crossa, 1990). The AMMI model combines the analysis of variance for the genotype and environment main effects with principal component analysis of the genotype x environment interaction (Zobel et al., 1988). The AMMI analysis is hence most widely used method for GEI investigation as it captures a large portion of the G x E sum of squares; it clearly separates main and interaction effects and often provides meaningful interpretation of data to support breeding programmes such as genotype stability (Gauch and Zobel, 1996). 
       
Each stability parameter has its strengths and weaknesses for the selection of desirable varieties and specific ways for addressing GEI effects; therefore it was recommended by several workers to incorporate both parametric and non-parametric methods (Becker and Leon, 1988; Mohammadi and Amri, 2008). To identify high-yielding and stable genotypes average of sum ranks (ASR) of parametric and non-parametric measures can be used as an effective strategy (Vaezi et al., 2019). Therefore, in the present study an effort had been made to identify stable and high yielding pigeonpea genotypes for variable environmental conditions by integrating various parametric and non-parametric models of stability.

Plant materials and field evaluation
 
Field trials consisting of twenty pigeonpea genotypes including 18 newly developed high yielding genotypes and two released varieties (Table 1) were laid down in Randomized Block Design with three replications at the N. E. B. Crop Research Centre, GBPUA&T, Pantnagar, Uttarakhand for four consecutive years in same plot. The years were treated as environments [Environment I (2016-17), Environment II (2017-18), Environment III (2018-19) and Environment IV (2019-20). The experimental site geographically falls in the humid sub-tropical climate zone situated at 29.5° North latitude, 79.3° East longitudes at an altitude of 243.84 m above mean sea level. The climatic condition remains highly variable in this region during different years (Fig 1). All the recommended package of practices was adapted to raise a normal and healthy crop. At harvest, seed yield data was collected for each plot and converted to kg/ha.
 

 

 
Statistical analysis and procedures
 
An online program, STABILITYSOFT (Pour-Aboughadareh et al., 2019) was used to calculate the parametric and non-parametric stability measures.  Parametric measure includes Wricke’s (1962) ecovalence (W²_i); Shukla’s (1972) stability variance (σ_i²); and Francis and Kannenberg’s (1978) coefficient of variance (CVi). Non-parametric measures include Huehn’s (1990) and Nassar and Huehn’s (1987) non-parametric measures (S⁽ⁱ⁾); Thennarasu’s non-parametric (1995) measures (NP⁽ⁱ⁾) and Kang’s (1988) rank-sum. The regression coefficient (b_i) and deviation from regression (S²d_i) parameters (Eberhart and Russell, 1966) were estimated by using web based statistical tool PBSTAT-GE 2.3 (www.pbstat.com). All the statistical analysis related to AMMI model as well as constructions of biplots were performed on GEA-R (2017) Version 4.1 software available at www.cimmyt. org. AMMI stability values (ASV) were also estimated to rank the genotypes according to their stability (Purchase et al., 2000). The Yield Stability Index incorporating both mean yield and stability in a single index was used to identify stable genotypes with high seed yield (Kang, 1993; Bajpai and Prabhakaran, 2000; Bose et al., 2014). The AMMI I, AMMI II biplots were obtained by using the software PBSTAT-GE 2.3. A hierarchical cluster analysis (HCA) was performed to group the tested genotypes. The square Euclidean distance was used as the dissimilarity measure required in Ward’s clustering method (Ward, 1963). The statistical analyses were performed using SPSS (2007) version 16.0 software packages.
 

The data obtained from studied environments were analysed separately to check the significance of different genotypes. The environments where significant genotypic differences were observed were used further for pooled analysis. A close perusal of Table 2 indicated that in the present study genotypic differences were significant under all the studied environments and hence all the four environments were used for pooled analysis. The pooled analysis was done by using the mean data of four environments and it was found that the mean sum of squares for genotype, environment and G x E interaction were highly significant (p<0.01) (Table 3). The presence of significant genotypic differences for seed yield among the genotypes under the studied environments indicated the preponderance of sufficient genetic variability among the genotypes. The mean seed yield over all the environments ranged from 853.27 kg/ha (PA 617) to 1774.85 (PA 622) with a grand mean of 1231.456 kg/ha (Table 1). The significant differences among different environments indicated that these environments are different in their conditions and prevailing climatic conditions influenced the seed yield to a large extent. The significant G x E interaction indicated that genotypes performed differently under different environments. The significance of genotype, environment and G x E interaction effects for seed yield in pigeonpea was reported earlier by several researchers (Meena et al., 2017, Muniswamy et al., 2018; Singh et al., 2018 and Gaur et al., 2020). As the G x E interaction was found significant, the analysis was proceeded further to estimate stability parameters by different models.
 

 

 
Stability analysis by using parametric models
 
Wricke’s (1962) ecovalence model (W²_i) evaluates stability on the basis of the contribution of each genotype to the total G x E sum of squares. The ecovalence (W²_i) of a genotype is explained as its interaction with the environments, squared and summed across environments. The genotypes having a low value of W²_i are considered as stable. The study of ecovalence (W²_i) indicated that genotype UPAS 120 (W²_i= 11647.56, rank=1) was most stable followed by genotype PA 620 (W²_i=12278.52, rank=2) and PA 622 (W2i=18185.33, rank=3) whereas, the genotype PA 631 (W²_i=345050.60, rank=20) was found as least stable (Table 4).  Shukla (1972) suggested that variance component of each genotype across environments can be a reliable measure of stability. A stable genotype exhibited a stability variance (σ²) either equal to zero or a value non-significantly deviates from zero. The obtained values of stability variance (σ²) indicated that UPAS 120 (σ²= 2490.77) was most stable genotype followed by PA 620 (σ²= 2724.46) and PA 622 (σ²= 4912.17) whereas the genotype PA 631 was least stable (σ²=125973.40).  The Table 4 indicated that the ranking of genotypes by Shukla’s stability variance (σ²) is same as Wricke’s ecovalence (W²_i). This similarity in ranking by these two methods is due to the fact that the stability variance (σ²) is a linear combination of the ecovalence (W²_i) and therefore both W²_i and σ² are equivalent for ranking purposes. As a result, it is adequate and acceptable to use either of these two statistics.  In case of coefficient of variation (CV_i) proposed by Francis and Kannenberg (1978) the genotypes with low CV_i and high mean yield are considered to be the most desirable. A close perusals of Table 4 indicated that genotype PA 620 (CV_i= 6.657, rank=1) followed by PA 629 (CV_i= 7.052, rank=2) and PA 622 (CV_i= 8.678, rank=3) were most stable genotypes whereas the genotype PA 621 (CV_i= 26.044, rank=20) was found to be least stable. The rank of the genotypes by parameter CV_i is different as compared to the W²_i and σ² stability parameters. The major limitation of this model as pointed out by Bowman and Watson (1997) is that  in comparing genotypes across high and low yielding environments, if the mean and standard deviation do not vary in a parallel way as performance increases, a bias would happen, whereby high mean values result in low CV and low mean values in high CVs.  The model proposed by Eberhart and Russell (1966) was one of the most preferred methods of stability analysis because of its simplicity and reliability. In Eberhart and Russell (1966) model, the stable genotypes were identified based on high mean yield, regression coefficient (b_i) around unity and mean square deviations from regression (s²d_i) non-significant from zero. The results indicated that only one genotype i.e. PA 620 (b_i=0.98, rank=1; s²d_i =3928.65, rank=1) performs better across all studied environments and hence considered as most stable.  The similar kind of results by using Eberhart and Russell (1966) model in different pigeonpea genotypes  were also obtained by Reddy et al., (2011) and Patel and Tikka (2014).
 

 
Stability analysis on basis of AMMI biplots, ASV and YSI
 
The ANOVA of AMMI revealed that for grain yield, the environment, genotype and G x E interaction was found to be significant (Table 5). This indicated that seed yield was influenced by both, main effects as well as their interactions.  In pigeonpea, the significance of main effects (environment and genotype) as well as G x E interaction effects for seed yield was also reported earlier by (Chauhan et al., 1999; Wamatu and Thomas, 2002; Muniswamy et al., 2018; Singh et al., 2018 and Gaur et al., 2020). An insight of Table 6 indicated that for seed yield per plant, 72.12 % of total sum of square (TSS) was attributable to genotypic effects, 22.57 % to G x E effects and 5.30 % to environment effect. The IPCA I, IPCA II and IPCA III accounted for 60.27 %, 28.02 % and 11.69 % of the genotype x environment interaction respectively.  The major portion of total sum of squares (TSS) was contributed by genotypic effects indicating the preponderance of genetic diversity in the genotypes under study. The major contribution of genotypes towards the total sum of squares for seed yield per plant and other important components has also been reported earlier by Singh et al., 2018; Sharma et al., 2019 and Gaur et al., 2020. The effects by environment were small as compared to genotype and G x E interaction effects but still exhibited significance indicating that the environments under study were variable. The larger contribution of G x E interaction than the environment suggested the differential response of environments towards genotypes.  The sum of squares due to G x E interaction were further partitioned into three principal component axis i.e. IPCA I, IPCA II and IPCA III accounting for 100 per cent of the G x E interaction sum of squares and uses entire degrees of freedom available in the interaction. Thus, in the present study, AMMI having three principle components axis was found as the best predictive model.
 

       
AMMI biplots provide a visual inspection and interpretation of the G x E interaction (Gabriel, 1971). In AMMI I biplot IPCA 1 scores close to zero have small interactions and therefore, considered more stable (Vargas and Crossa, 2000). In AMMI biplot II, the genotypes near the origin are non-sensitive to environmental interaction, hence are more stable and those distant from origin are sensitive and have large interaction. On the basis of AMMI model (IPCA I and IPCA II components) the AMMI stability value (ASV) can be calculated (Purchase et al., 2000). On the basis of AMMI biplot I and II, ASV (AMMI Stability Value), genotype PA 620 (IPCA I, 0.02912; ASV rank=1) was identified as most stable genotype followed by UPAS 120 (IPCA I, -0.04768; ASV rank=2) and PA 622 (IPCA I, -0.14601; ASV rank=3) (Table 6 and 4 and Fig 2-3).  Stability per se should not be used as the sole selection criteria to identify desirable genotypes as most of stable genotypes would not necessary give the high seed yield (Mohammadi et al., 2007). The Yield Stability Index (YSI) is an integrated approach that incorporates both mean performance as well as stability in a single index and can be used for simultaneous selection of high yield and stability (Kang, 1993; Bajpai and Prabhakaran, 2000). On basis of Yield Stability Index (YSI) scores, the genotype PA 620 (YSI=4, rank=1) followed by PA 622 (YSI=4, rank=1) and PA 626 (YSI=9, rank=2) were identified as most stable and high yielding genotypes for seed yield per plant across studied environments (Table 5 and 4).
 

 

 

 
Stability analysis by using Non-parametric models
 
Huehn (1990) and Nassar and Huehn (1987) suggested four non-parametric statistics i.e. S⁽¹⁾, the mean of the absolute rank differences of a genotype over all tested environments, S⁽²⁾, the variance among the ranks over all tested environments, S⁽³⁾, the sum of the absolute deviations for each genotype relative to the mean of ranks and S⁽⁶⁾, the sum of squares of rank for each genotype relative to the mean of ranks. The lowest value for each of these statistics reveals high stability of genotype.  The perusals of Table 7 indicated that the stability parameter S⁽¹⁾ clearly distinguished stable and unstable genotypes. Genotype PA 622 (S⁽¹⁾ = 0.5, rank=1), followed by PA 620 (S⁽¹⁾ = 1.16, rank=2) and PA 617 (S⁽¹⁾ = 1.5, rank=3) were found as most stable genotypes whereas the genotype PA 631 (S⁽¹⁾ = 8, rank=20) was least stable.  The stability parameter S⁽²⁾ also exhibited same result as obtained by S⁽¹⁾. The genotype PA 622 (S⁽²⁾ = 0.25, rank=1), followed by PA 620 (S⁽²⁾ =0.91, rank=2) and PA 617 (S⁽²⁾ =1.58, rank=3) were found as most stable genotypes whereas the PA 631 (S⁽²⁾ = 40, rank=20) was least stable. The parameter S⁽³⁾ also exhibited similar kind of results like  S⁽¹⁾ and S⁽²⁾ parameters the most stable genotype was identified as PA 622 (S⁽³⁾ = 0.037, rank=1) followed by PA 620 (S⁽³⁾ =0.154, rank=2) but here UPAS 120 (S⁽³⁾ = 0.574, rank=3)  is identified as third most stable genotype in place of PA 617 and the genotype PA 621 S⁽³⁾ =10.444, rank=20) was found as least stable in place of genotype PA 631.  In case of stability parameters S⁽⁶⁾ the genotypes PA 622 (S⁽⁶⁾ = 0.075, rank=1), followed by PA 620 (S⁽⁶⁾ = 0.169, rank=2) and PA 619 (S⁽⁶⁾ = 0.318, rank 3) were found as most stable. The use of S^(1-6) models suggested that PA 622 is the most stable genotype as it ranked one in all four S⁽ⁱ⁾ parameters followed by PA 620 which ranked second in all four S⁽ⁱ⁾ parameters indicating the equal efficiency of these models and any one of these parameters can be used. Four NP^(1–4) statistics are a set of alternative non-parametric stability statistics defined by Thennarasu (1995). These parameters are based on the ranks of adjusted means of the genotypes in each environment. Low values of these statistics reflect high stability. The stability parameter NP⁽¹⁾ revealed UPAS 120 (NP⁽¹⁾ = 2, rank=1) and PA 622 (NP⁽¹⁾ = 2, rank=1) followed by PA 620 (NP⁽¹⁾ =2.25, rank=3) as most stable whereas the genotype PA 618 (NP⁽¹⁾ = 8, rank 20) was least stable. The stability parameter NP⁽²⁾ revealed PA 628 (NP⁽²⁾ = 0.153, rank=1), PA 624 (NP⁽²⁾= 0.160, rank =2 ) and PA 616 (NP⁽²⁾=0.166, rank=3) as most stable while the genotype PA 617 (NP⁽²⁾=2.916, rank=20) was least stable. The stability parameter NP⁽³⁾ revealed PA 620 (NP ⁽³⁾ = 0.134, rank=1), PA 622 (^NP(3)= 0.136, rank=2) and UPAS 120 (NP⁽³⁾= 0.229, rank=3) as most stable while the genotype  PA 617 (NP⁽³⁾= 2.040, rank=20) was least stable. The stability parameter NP ⁽⁴⁾ revealed PA  622 (NP⁽⁴⁾= 0.025, rank=1) PA 620  (NP⁽⁴⁾= 0.065, rank=2 ) and PA 619 (NP⁽⁴⁾= 0.125, rank=3) as most stable while the genotype PA 632  (NP⁽⁴⁾= 0.8, rank=20) was least stable. The NP⁽¹⁾ and NP⁽⁴⁾ revealed PA  622 as the most stable genotype while NP ⁽³⁾ gave it second rank and NP⁽²⁾ gave it eight rank. These results suggested that NP⁽²⁾ is least efficient among all the NP⁽ⁱ⁾ statistics while NP⁽¹⁾ and NP⁽⁴⁾ are most effective. Kang’s rank-sum (Kang, 1988) used both yield and stability statistics to identify high-yielding and stable genotypes. The genotype with the highest yield and lower σ²_i are assigned a rank of one. Then, the ranks of yield and stability variance are added for each genotype and the genotypes with the lowest rank-sum are the most desirable.  The analysis of stability by using Kang rank sum reveals that the genotype PA 622 is most stable as it has Kang rank 1 stable followed by PA 620 ( rank=2) and UPAS 120 (rank=3) whereas, the genotypes PA 615 and PA 621 were least stable as both of them has rank of 19.
 
@table 7
       
The results obtained revealed that separate application of parametric and non-parametric models results in differential ranking of genotypes which creates an ambiguity in selection of high yielding and stable genotype. As in case of parametric model, most of the parameter suggested that genotype PA 620 was the most stable genotype (CV_i, b_i, s²d_i, ASV and YSI), however, the W²_i and σ² parameters suggested that genotype UPAS 120 was most stable genotype.  Similarly in case of non-parametric model the genotype PA 622 was found as most stable (S^(1-6), NP⁽¹⁾, NP⁽⁴⁾  and KR) however the parameter NP⁽³⁾ gave it second rank and NP⁽²⁾ gave it eight rank. Due to the difficulty in identifying high-yielding stable genotypes based on a single stability measure, the use of average of sum of ranks (ASR) of all measures to select high yielding and stable genotypes with low ASR values was recommended (Vaezi et al. 2019). The perusal of Table 5 indicated that the genotypes PA 622 and PA 620 were most stable and high yielding genotypes as they had lowest ASR value of 2.0 and 2.18 respectively. UPAS 120 was also showing good stability as it had an ASR value of 2.87. The genotype PA 621 was the least stable genotype as it had highest ASR value 15.75.
               
In order to confirm the above findings, HCA based on seed yield and ASR ranks was used and it grouped the genotypes into two main clusters i.e. Cluster I (CI) and Cluster II (CII) (Fig 4). Cluster I comprised of low yielding genotypes with relatively high values of ASR i.e. PA 630 (yield=987.50 kg/ha, ASR=13.43), PA 632 (yield=971.22 kg/ha, ASR=13.43), PA 615 (yield=891.82 kg/ha, ASR=12.75), PA 627 (yield=991.67 kg/ha, ASR=12.62), PA 621 (yield=1153.15 kg/ha, ASR=15.75), PA 617 (yield=853.27, ASR=10.68), PA 623 (yield=890.82 kg/ha, ASR=9.81) and PA  629 (yield=874.92 kg/ha, ASR=7.87), suggesting that these genotypes may have specific adaptations to some of the environments. Cluster II was divided into two sub-clusters (sub-cluster III and sub-cluster IV), sub-cluster III mainly comprised high to moderate-yielding genotypes with low ranks of stability [PA 620 (yield=1579.92 kg/ha, ASR=2.18), PA 622 (yield=1774.85 kg/ha, ASR=2.20), PUSA 992 (yield=1331.17 kg/ha, ASR=5.68), PA 626 (yield=1571.40 kg/ha, ASR=5.56), PA 628 (yield=1271.50 kg/ha, ASR=7.06), UPAS 120 (yield=1268.57 kg/ha, ASR=2.87)] indicating that these genotypes were high yielding along with high stability across the environments. Hence, these genotypes could be effectively used to improve performance and adaptation in pigeonpea breeding programs. The sub-cluster IV consisted of high to moderate yielding genotypes with high ASR values [PA 616 (yield=1232.57 kg/ha, ASR=8.06), PA 624 (yield=1259.47 kg/ha, ASR=8.37), PA 625 (yield=1232.85 kg/ha, ASR=11.25), PA 618 (yield=1552.40 kg/ha, ASR=11.56), PA 631 (yield=1330.17 kg/ha, ASR=14.81), PA 619 (yield=1609.82 kg/ha, ASR=8.68)] suggesting that these genotypes performs well only under certain environments.

From the above findings it can be concluded that parametric and non-parametric models results in differential ranking of genotypes. Genotype PA 620 ranked first in most of the parametric models and genotype PA 622 rank first in case of most of the non-parametric model. Hence it is better to integrate these models and use the ASR method for a more precise selection of high yielding and stable genotypes. The genotypes having low ASR values can be recommended for cultivation across a range of environments. The use of HCA based on mean yield and ASR method indicated that the genotype PA 622 (yield=1774.85 kg/ha, ASR=2.00), PA 620 (yield=1579.92 kg/ha, ASR=2.18), UPAS 120 (yield=1268.57 kg/ha, ASR=2.87), PA 626 (yield=1571.40 kg/ha, ASR=5.56), PUSA 992 (yield=1331.17 kg/ha, ASR=5.68) and PA 628 (yield=1271.50 kg/ha, ASR=7.06) as most stable and high yielding.