8. Case study: durum wheat in Algeria

8.1 Experiment data

Layout of experiments

Within the framework of the bilateral cooperation project “Amélioration et renforcement du système nationale d’adaption variétale du blé dur en Algérie” funded by the governments of Algeria and Italy and carried out by the Institut Technique des Grandes Cultures (Algeria) and the Istituto Agronomico per l’Oltremare (Italy), 24 durum wheat cultivars of different origin were grown for two years on 18 Algerian locations. They comprised: 19 improved varieties (of which 10 were released from international research centres, 3 from Algeria, 3 from Italy, and 1 each from France, Spain and Tunisia); 4 Algerian landraces; and 1 old cultivar of French origin. The trials aimed to:

• optimize the variety recommendation across the country for this fundamental staple crop; and

• support the decisions on adaptation and yield stability targets of the national breeding programme.

Detailed information on experiments and major results of the project are reported by Annicchiarico et al. (2002a, 2002b, 2002c). Each trial was designed as a randomized complete block with four replicates. One site suffered two crop failures and was, therefore, excluded from the analysis. The geographic position of the remaining sites is shown in Figure 8.1, while altitude and mean yield values are reported in Table 8.1. Site ‘13’ experienced a crop failure in the first year; second year data for sites ‘3’ and ‘15’ were not available. Therefore, initial ANOVAs and analysis of adaptation of grain yield values concerned the balanced data set composed of 14 locations by 2 years. Information from sites with one year’s data was exploited at a later stage of the analysis.

FIGURE 8.1 - Geographic position of test locations and definition of two subregions for a specific adaptation strategy based on pattern analysis and AMMI + cluster analysis results

Source: Annicchiarico, 2002c.

Data files

The following IRRISTAT data files were progressively created and used for statistical analyses:

i) Data file A (for ANOVAs), including original and log₁₀-transformed plot values of the 14 test sites and holding genotype, location, year, block and environment (identifying single trials) as classification criteria.

ii) Data file B (for analysis of adaptation), including, for each genotype-location combination relative to the 14 sites: mean over years and replicates (obtained as an output file from IRRISTAT’s ANOVA module); genotype mean yield across locations, site mean yield, and site mean value across test years of each of six available environmental variables (assigned to each genotype-location combination through IRRISTAT’s logical operators); GL interaction effects (assigned to each genotype-location combination by the formula [4.1] through IRRISTAT’s logical operators); and a few empty variables for storing expected (fitted) values for factorial regression models.

iii) Data file C (for ANOVAs), similar to file A but including also plot values for the three sites with one year’s data.

iv) Data file D (for correlation and multiple regression analyses concerning site values), including mean yield, mean values of environmental variables, scaled PC scores for original and transformed data, a few dummy variables used for discriminant analysis, and a few empty variables for storing expected values for regression models, for all test sites.

v) Data file E (for location classification by cluster analysis), including the following four variables: location, number of the variable for site classification (with levels 1 and 2), scaled PC scores of sites for transformed data and one environmental variable of sites (since AMMI-1 and one-covariate factorial regression models were selected, values for the second classification variable relative to PC 2 or a second covariate corresponded to a dummy variable with extremely limited range of variation).

vi) Data file F (for genetic correlation analysis), including original yield values of genotypes on the three sites with one year’s data, and mean values of genotypes across sites for each of two subregions.

vii) Data file G (for worksheet calculation of expected AMMI- and factorial regression-based nominal yields or nominal yield reliabilities as a function of fixed levels of two climatic variables), including rainfall and winter mean temperature as variables and all possible combinations of 10 rainfall by 10 temperature levels as observations (rows).

In addition, one Excel data file (obtained by exporting IRRISTAT’s Data file B) aided the calculation of within-location phenotypic variance of genotype values (the calculation of which could easily be performed by IRRISTAT using the Data selection option now implemented in the Summary module).

TABLE 8.1 - Code, altitude, mean yield, scaled score on the first GL interaction PC axis for original and log₁₀-transformed yield data, winter mean temperature and rainfall amount of test sites

Note: See Figure 1 for geographic position of sites.

^a Values into brackets relate to only one test year.
^b Average daily temperature across January and February months of test years.
^c From September to June.

Source: Annicchiarico et al., 2002a.

Equations relative to scaling-up of results were calculated using the software adopted for GIS management, to obtain the desired graphical output. All remaining statistical analysis or calculation, with the exception of a limited amount of manual calculation, was performed by IRRISTAT (to which, therefore, all cited analytical modules apply), in accordance with previous indications provided in Sections 4.5, 5.9 and 7.3. Unless specified, GL interaction PC scores are hereafter intended as scaled scores (i.e. u’ and v’ in previous notations).

8.2 Adaptation strategy and yield stability targets for breeding

Lack of genetic correlation and heterogeneity of genotypic variance

Analysis for this objective follows the flow chart in Figure 2.3 and the relative discussion in Section 2.6. The performance of:

• one ANOVA per individual environment (applying the Single Site Analysis module with the Multiple Selection option on Data file A); and

• a combined ANOVA including the random factors, genotype and environment, in addition to block within environments (applying the ANOVA module on Data file A),

and the application of formulae reported in Section 4.3 and Table 4.1, allowed for estimation of the variance components relative to genotype, GE interaction, heterogeneity of genotypic variance among environments, and lack of genetic correlation among environments. The last component (0.047 [t/ha]²) proved about as large as the genotypic component (0.046 [t/ha]2), justifying further analysis (Fig. 2.3). The estimated pooled genetic correlation among environments (r_g = 0.49) confirmed that GE interaction effects are considerable. However, the larger size of the heterogeneity of genotypic variance (0.063 [t/ha]²) relative to the lack of genetic correlation component pointed to the transformation of data.

The criterion described in Section 5.6, based on the regression of the within-site phenotypic variance of genotype values on the site mean yield expressed on a logarithmic scale (applying the Regression module), supported (for b greater than zero [P < 0.01] and equal to 1.92) the adoption of a logarithmic transformation for most subsequent analysis (performed on kg/ha rather than t/ha yield units to avoid negative transformed values). In fact, the estimation of the GE interaction components of variance for transformed data (not strictly needed) showed a substantial reduction in the heterogeneity of genotypic variance relative to lack of genetic correlation component (0.0021 relative to 0.0041 [log₁₀ kg/ha]²).

Variance of genotype, GL interaction and other GE interaction effects

A combined ANOVA performed on transformed data and including, besides the block factor, genotype, location and year as crossed factors (applying the ANOVA module on Data file A), and the formulae reported in Table 4.2 for Model 1 of ANOVA, allowed for the estimation and testing of the variance components relative to genotype and GL, GY and GLY interactions. DF and MS values for all sources of variation are reported in Table 8.2. An appropriate error MS for the genotype factor and the associated DF were estimated from ANOVA information through the formulae reported in Section 4.1, obtaining, for transformed data, the following values: M_err = 0.0389; DF ≈ 28. The genotypic, GL interaction and GLY interaction components of variance were significantly different from zero, unlike the GY interaction variance (Table 8.2). The estimated variance of the GL interaction was large enough relative to the genotypic variance (≈ 35%) to justify the analysis of adaptation (see Section 2.6).

Analysis of adaptation

Table 8.2 also reports joint regression and AMMI analysis results obtained by analysing genotype by location cell means of transformed data (in Data file C) by the G × E Interaction module. For integration with ANOVA results, SS values outputted by IRRISTAT for individual PC axes, heterogeneity of genotype regressions and residual GL interaction terms were multiplied by N’ = no. years × no. replicates = 8, prior to calculation of MS values. The application of the F_R test to PC axes implied the same P level as the GL interaction term for PC 1, and the following calculation for PC 2 (see Section 5.3):

For PC 1:

DF = 24 + 14 - 1 - 2 = 35

For the residual GL interaction from AMMI-1:

DF = 299 - 35 = 264
SS = 10.375 - 5.407 = 4.968
MS = 4.968/264 = 0.0188
F ratio = 0.0188/0.0285 = 0.66

Since the F ratio was lower than the critical value F ≈ 1.25 (relative to P < 0.05 and DF = 264, 299), the selected AMMI model included only PC 1.

Despite the significance of heterogeneity of regressions and the lack of significance of deviations from regressions, the joint regression model showed a rather modest ability to explain GL effects both in absolute terms (R² = 21%) and in comparison with AMMI-1 (R² = 52%). In this and subsequent analyses, the sign of site and genotype PC 1 scores as outputted by the program was consistently reversed to facilitate the comparison with joint regression modelling, thereby assigning positive and negative score values to locations tending towards high and low mean yield, respectively (the double change of sign has no implication for the results - see Section 5.3). However, site score on PC 1 was only moderately correlated with site mean yield expressed either in original or transformed units (r < 0.60, P < 0.05), while it was strictly correlated (P < 0.001) with two environmental variables, namely, mean temperature in winter (r = 0.86) and altitude (r = -0.83) of sites (for analyses performed by the Regression module on Data file D). Site values for these variables are reported in Table 8.1.

TABLE 8.2 - ANOVA results and estimate of variance components, for original and transformed yield of 24 genotypes grown at 14 locations for two years

Note: GL interaction partitioned by joint regression (JR), AMMI and factorial regression (FR) models.

^a ns = not significant; *, ** and *** = significant at P < 0.05, P < 0.01 and P < 0.001, respectively.

^b Heterogeneity of genotype regressions summarizes 45% and 21% of the GL interaction SS for original and transformed data, respectively.

^c The first PC axis summarizes 70% and 52% of the GL interaction SS for original and transformed data, respectively.

^d Best models among those including up to five covariates. Heterogeneity of genotype regressions summarizes 61% and 40% of the GL interaction SS for original and transformed data, respectively.

Source: Annicchiarico et al., 2002a, 2002c.

Factorial regression results are also reported in Table 8.2. One-covariate models relative to each of the six environmental variables were assessed by analysing the GL interaction effects for transformed data (in Data file B) as a function of each covariate (indicated as the Site Index in the G × E Interaction module). The model with the highest R² value included winter mean temperature (R² = 40%). For insertion in Table 8.2, the outputted SS for the covariate was multiplied by N’ = 8. The derived MS was tested on the GLY interaction MS.

For each possible two-covariate model including the best single covariate, multiple regression analyses were performed separately for each genotype on the GL interaction effects for transformed data (applying the Regression module with the Multiple data selection option on Data file B). Both the regression SS and the deviation from regression SS values were summed up over analyses and multiplied by N’ = 8 to obtain, respectively, the model SS and the residual GL interaction SS. No additional covariate was significant after calculating and testing its MS. For example, the MS value for rainfall amount was calculated from its partial regression SS, equal to the difference between the two-covariate SS (4.6120) and the one-covariate SS relative to winter mean temperature (4.1288), divided by DF = 23.

The variance of significant GL interaction components (estimated as described in Section 5.4 for comparing analysis of adaptation models) is reported in Table 8.2. For example, the calculation for heterogeneity of genotype regressions in the joint regression model was:

S_C² = (0.0953 - 0.0285)/(4 × 2 × 14) = 0.0006 (log₁₀ kg/ha)²

This criterion suggested an advantage of AMMI-1 and factorial regression models over joint regression. Also the criterion proposed by Brancourt-Hulmel et al. (1997) provided the same indication - the ratio of % GL interaction SS to % GL interaction DF being 2.74, 4.45 and 5.17 for joint regression, AMMI-1 and factorial regression, respectively. For example, the ratio value for the first model derived from the fact that it accounted for 21.1 and 7.7 percent of the GL interaction SS and DF, respectively. As a result of the comparison, the joint regression approach may be neglected in further analyses.

Provisional definition of subregions

The next step of analyses (Fig. 2.3) is exemplified for all methods discussed in Sections 5.2 to 5.5. Test sites were partitioned into two groups on the basis of the estimated main crossover point for each model on the one hand, and by cluster analysis on the other. Formulae [5.6] and [5.8] for estimating main crossover points required information on genotype parameters that was either outputted from previous applications of the G × E Interaction module (site scores on PC 1, and factorial regression b_i values of genotypes), or calculated from other information also outputted by the same applications (b_i values in formula [5.6], and intercept values a_i in both formulae - see Section 8.3 for details of nominal yield calculation). The position of the main crossover point in the site ordination for PC 1 based on transformed data is reported in Figure 5.9.

Cluster analysis using Ward’s method was performed on site PC 1 scores on the one hand, and site mean value of winter mean temperature on the other (applying the Pattern analysis module on Data file E, for classification based on raw data of either variable plus a dummy, uninfluent variable). Scaled PC 1 scores might have been converted to original PC 1 scores through multiplication by the square root of the singular value for PC 1 outputted by the program (l₁ = 0.8220; √ l₁ = 0.9067), but the transformation has no influence on results when only one PC axis is involved in the cluster analysis. The classification based on the site covariate could be extended to the three test sites with missing data in one test year (since covariate values were available for both test years). The results, using the lack of significant GL interaction within groups of locations (P < 0.05) as the truncation criterion (verified for sites with two years’ data by separate ANOVAs, applying the ANOVA module with the Data selection option on Data file A), supported the definition of two groups of locations for both cluster analyses.

Site classification was also carried out by pattern analysis, applying the Pattern analysis module on genotype by location cell mean values of original yield (in Data file B) with the default options (which imply the within-location standardization of genotype values, and the adoption of Ward’s method for cluster analysis). The above truncation criterion, requiring separate ANOVAs for groups of locations and performed, as earlier, on log-transformed data, supported also in this case the definition of two subregions (in fact, the ANOVAs performed on original data would rather have supported the definition of three subregions).

Specific adaptation advantage

Once the provisional subregions were defined, wide vs. specific adaptation strategies were compared (Fig. 2.3) on the basis of yield gains predicted from the same data set (as described in Section 6.1). For this case study, results are also reported for other classification criteria, which appeared less useful in the analyses and would not, therefore, be considered for provisional definition of subregions, such as those based on joint regression analysis or analysis of adaptation performed on original yield data (Table 8.3; no reference is made, however, to data sets and IRRISTAT procedures needed for these additional analyses). Since the factorial regression modelling of original data identified two significant environmental variables (Table 8.2), the cluster analysis procedure that complemented it could also contemplate the assignment to each variable of a weight that is proportional to its partial regression SS in the regression model, whereas site classification by the main crossover criterion was not applicable. The scenario with two subregions was retained in all cases, to facilitate the comparison of classification criteria.

Prior to the assessment, the three test locations with missing data in one test year were allocated to either subregion by different procedures depending on the classification criterion (see Section 5.7). For criteria based on AMMI or pattern analysis, each site was assigned to the subregion with which it showed the highest genetic correlation. The correlation was estimated from analyses on original yield values by the procedure described in Section 5.7, which implied:

• an ANOVA for each of the three test sites (applying the Single site analysis module with the Data selection option on Data file C);

• a separate ANOVA, including the factors genotype and environment besides block, for each of the two groups of sites (applying the ANOVA module with the Data selection option on Data file A); and

• the estimation of phenotypic correlations of genotype values across each subregion with those in each of the three test sites (applying the Regression module on Data file F), besides some manual calculation for estimation of variance components and broad sense heritability values.

TABLE 8.3 - Comparison of analytical methods for definition of two subregions for durum wheat breeding in Algeria. Predicted yield gain over the region from a specific adaptation strategy (DGS) relative to wide adaptation (DGW)

^a JR = joint regression; AMMI = model with one GL interaction PC axis; FR = factorial regression holding winter mean temperature alone (transformed data) or with annual rainfall (original data) as covariate; CA = cluster analysis (-nw = not weighted, and -w = weighted, covariates in the FR-based approach); CP = main crossover point.

^b P = proportion of the target region; E = number of selection environments; ΔG = predicted yield gain per selection cycle (four experiment replicates; standardized selection differential = 1.64, for 10% selection intensity applied to 20 elite breeding lines).

^c ΔG_W = 0.297 t/ha, from selection in E = E_A + E_B = 6 environments.

Source: Annicchiarico et al., 2002c.

Of the 13 classification criteria (Table 8.3), only those two that were based on joint regression analysis of transformed data provided an identical classification of test locations (data not reported).

The two adaptation strategies hypothesized in all cases:

• three sites and two years (i.e. six environments) for final stages of selection, of which two sites assigned to the larger subregion and one to the smaller in the specific adaptation prospect (as exemplified in Fig. 6.1);

• four replicates per selection environment; and

• 10% selection intensity applied to 20 breeding lines, i.e. standardized selection differential = 1.64 (Falconer, 1989^[35]).

The average yield gain across the region (ΔG_W) was predicted by formula [6.2], after estimating the relevant components of variance through a combined ANOVA performed on original yields of all (i.e. 31) test environments (applying the ANOVA module on Data file C) and the application of formulae for Model 1 in Table 4.1. In particular, from the estimates s_g² = 0.0470, s_ge² = 0.1039 and s_e² = 0.0915 (t/ha)², the estimates h_AB² = 0.690 and s_p(AB) = 0.262 were obtained using formula [6.1] and the square root of its denominator, respectively (placing E = 6), with an estimated gain ΔG_W = 0.297 t/ha per selection cycle.

As an example of predicted yield gain derived from specific breeding, consider the two subregions A and B provisionally identified by AMMI + cluster analysis of transformed data. Subregion A (including seven out of 17 test sites) was assigned the proportion of the target region P_A = 0.41 and one selection site (i.e. E_A = 2), whereas subregion B was assigned P_B = 0.59 and two selection sites (E_B = 4). Components of variance relative to genotype, GE interaction and pooled error were estimated for each subregion through formulae for Model 1 in Table 4.1, after performing a separate ANOVA for each of the two groups of test environments (applying the ANOVA module with the Data selection option on Data file C). Then, heritability and phenotypic standard deviation values were estimated with formula [6.1] (with E = 2 or E = 4 as appropriate) and inserted in formulae for prediction of yield gains for subregions A and B. Finally, these yield gains were used for predicting the average yield gain across the region (ΔG_S). In particular:

ΔG_A = i h_A² s_p(A) = 1.64 × 0.38 × 0.19 = 0.118 t/ha per selection cycle

ΔGB = i h_B² s_p(B) = 1.64 × 0.71 × 0.37 = 0.439 t/ha per cycle

ΔGS = (ΔG_A P_A) + (ΔG_B P_B) = (0.41 × 0.118) + (0.59 × 0.439) = 0.308 t/ha per cycle

Considering that a useful criterion for site classification tends to maximize the opportunities offered by a specific adaptation strategy, no criterion that would have been discarded in previous stages of the analyses ranked among the best ones according to the ΔG_S/ΔG_W ratio in Table 8.3. The possible definition of three subregions, suggested by some of these criteria, appeared even less useful for specific breeding (data not reported). On the whole, the results point to:

• the usefulness of the data transformation and the criterion underlying its adoption;

• no consistent advantage of cluster analysis over the main crossover criterion;

• no particular advantage of factorial regression-based criteria in this context; and

• the slight superiority of two methods, i.e. pattern analysis (using transformed data for the ANOVA-based truncation criterion) and AMMI + cluster analysis of transformed data.

Indeed, subregion definition using these latter two methods was very similar. With reference to Figure 8.1, the indications differed only for locations 7 and 14 (attributed to subregions A and B, respectively, by the former method, and to subregions B and A by the latter). Based on environmental data reported in Table 8.1, these sites represent a transitional zone separating the high-elevation, cold- and drought-prone subregion A from the warmer, partly drought-prone subregion B. While site mean yield is mainly associated with rainfall amount (Annicchiarico et al., 2002a), GL interaction effects relate mainly to winter cold, justifying the somewhat lower usefulness of the joint regression model.

A number of classification criteria have been considered in this comprehensive example. Just some of these criteria may be considered in ordinary applications, to reduce the calculation load. Pattern analysis, which is probably the least time-demanding criterion but one of the least informative on environmental and genotypic factors possibly contributing to GL interaction, could be recommended in all cases, in conjunction with a criterion based on AMMI or factorial regression analysis (depending on the available environmental data and the value of each analysis of adaptation model for the specific data set). A preference for joint regression-based criteria is justified only in those (probably infrequent) cases in which this model proves very accurate.

Scaling-up of results and adaptation strategy

The spatial and temporal scaling-up of results was attempted for subregion definition based on AMMI + cluster analysis and pattern analysis criteria, using GIS data. For example, analysis for the former criterion implied, with reference to formulae and notations in Section 5.8:

• the assignment of the dummy variable values Y_A = 0.5882 and Y_B = -0.4118 to sites of subregions A and B, respectively, calculated from the number of sites N_A = 7 and N_B = 10 for the two subregions;

• the performance of a stepwise multiple regression analysis of the dummy variable as a function of the environmental variables using a forward selection method (applying the Regression module with the Selection of variable option on Data file D), obtaining the following best model:

  Y_j = 1.2568 - 0.1528 V_j

  where V_j = mean winter temperature (as recorded in Table 8.1), with R² = 48%;

• the calculation of the dividing point for discrimination C = 0.0423 from the mean values for the two subregions, Z_A = 0.2819 and Z_B = -0.1973, of the expected Y_j values for the relevant test sites in test years (the Y_j values were outputted as fitted values by the Regression module);

• the calculation of Y_j values for new sites or test sites based on long-term GIS data of the environmental variable, and the allocation of sites to subregion A (if Y_j > C) and to subregion B (if Y_j < C).

For classification based on pattern analysis, the best model included only site altitude (R² = 46%). The results of discriminant analysis (very similar for the two classification criteria) provided a detailed description of subregion boundaries and modified the attribution to subregions of a few test sites characterized by unusually cold or warm winters in the test years (data not reported). However, the relative size of each subregion was modified. In particular, for the AMMI-based criterion the estimated relative size of the two subregions was: P_A = 0.345 and P_B = 0.655. Therefore, a better estimate of the predicted yield gain across the region provided by specific breeding (maintaining all previous hypotheses on selection environments, selection intensity etc.) was:

ΔG_S = (ΔG_A P_A) + (ΔG_B P_B) = (0.345 × 0.118) + (0.655 × 0.439) = 0.329 t/ha per cycle

which, compared with the gain predicted for a wide adaptation strategy (ΔG_W = 0.297 t/ha per cycle), indicated greater efficiency of specific breeding equal to: 0.329/0.297 = +10.6%. The advantage of this adaptation strategy was slightly lower (+8.7%) for subregion definition based on pattern analysis.

The advantage ultimately predicted for the specific adaptation strategy was probably underestimated by the fact that most tested cultivars are modern varieties bred at international research centres or in European countries and, as such, selected for wide adaptation or adaptation to areas which are quite different from the Algerian ones (especially the stressful subregion A). The subregion definition could conveniently be used for comparing adaptation strategies in terms of actual yield gains obtained from the selection of breeding lines across (wide prospect) or within (specific prospect) the two subregions.

Additional information

Further analysis (using procedures described in Sections 6.2 and 6.3) provided information on optimal selection sites, genetic resources and adaptive traits that could be exploited in breeding for wide or specific adaptation (Annicchiarico et al., 2002a, 2002b, 2002c). They highlighted, for example, the crucial importance of an appropriate matching of genotype phenology with the level and duration of cold stress in winter, and suggested that the two subregions require partly different plant ideotypes with respect to heading time, cold tolerance and plant stature (while being similar with respect to other useful traits). The information on selection sites could also be exploited to define a set of crucial test sites for the ordinary evaluation of novel varieties by public institutions.

Yield stability targets

The final step in the analysis regards the definition of yield stability targets. Should a wide adaptation strategy be preferred (implying parallel selection on a few sites representative of the two subregions), the choice between options 3 and 4 in Figure 2.3 depends on the size of the estimated GY and GLY components of variance for transformed data. Their summed value, equal to 239 percent of the genotypic variance (Table 8.2), supports the inclusion of yield stability as a breeding target (see Section 2.6).

If a specific adaptation strategy were preferred (as suggested by the results), decisions and relative analyses concern each subregion separately (option 5 vs. 6 in Fig. 2.3). Considering the subregion definition in Figure 8.1 (i.e. excluding from the assessment the two test sites placed in a transitional zone), the ANOVAs of transformed data performed for the two groups of environments, i.e. 11 environments for subregion A, and 16 for subregion B (applying the ANOVA module with the Data selection option on Data file C), and the application of formulae for Model 1 in Table 4.1, allowed for the estimation and testing of the genotypic and GE interaction components of variance. The latter - always different from zero (P < 0.001) - was distinctly larger than the former in subregion A (218%) but not in subregion B (110%). Therefore, yield stability may be considered a useful breeding target only for the former subregion, in which climatic conditions are not only more stressful but also more erratic, suggesting a higher number of selection environments or the choice of parent germplasm characterized by better temporal stability when breeding specifically for this subregion.

Correlation of yield stability values (expressed as the square root of environmental variance) with morphophysiological traits of material, for trait assessment limited to environments of subregion A, highlighted a few traits (e.g. taller plant stature and later heading time) associated with greater yield stability in this subregion.

8.3 Cultivar recommendation

Combined ANOVA

The analytical steps for the present objective follow the flow chart in Figure 2.4 and the discussion in Section 2.6. For reasons discussed earlier (Sections 4.4 and 5.6), all analyses were performed on original yield data.

TABLE 8.4. - Mean yield, environmental variance (S²) and Kataoka’s index of yield reliability (I) across environments, Type 4 stability variance (S_y(l)²), slope of regression on site mean yield (b), scaled score on the first GL interaction PC axis, and factorial regression equation for estimating GL effects, for eight genotypes

^a Data for 14 locations and two test years.

^b Data for 31 environments.

^c Mean separation of top-ranking values (letter ‘a’) according to Dunnett’s one-tailed test at P < 0.20.

^d *, ** = different from ‘GTA Dur’ at P < 0.05 and P < 0.01, respectively, according to Ekbohm’s test.

^e As lowest yield expected in 80% of cases.

^f *, ** = different from unity at P < 0.05 and P < 0.01, respectively.

^g Rainfall: from September to June (mm); temperature: as average of daily values in January and February (°C); *, ** = different from zero at P < 0.05 and P < 0.01, respectively.

Source: Annicchiarico et al., 2002a.

A combined ANOVA, with the factors: genotype, location, year and block, was performed according to Model 2 in Table 4.2 on the balanced data set relative to 14 locations (applying the ANOVA module on Data file A). Results are reported in Table 8.2. An appropriate error MS for the genotype factor and the associated DF were estimated through formulae reported in Section 4.1, obtaining the following values: M_err = 0.853; DF ≈ 37. From these values, Dunnett’s one-tailed critical difference for entry comparison over the target region may be calculated as d = 0.258 t/ha according to formula [4.2] (placing t’ ≈ 2.09 for DF ≈ 37, P < 0.20 and p = 24 entries in Table 4.6, and N = 112). Considering the subset of ten cultivars reported in Table 8.4, any entry that yields less than 1.94 t/ha is significantly lower yielding than the top-yielding genotype ‘GTA Dur’.

The occurrence of significant GL interaction in the ANOVA (Table 8.2) pointed, however, to the analysis of genotype responses to individual locations.

Analysis of adaptation and crossover interactions

The investigation of genotype responses is exemplified with respect to all analysis of adaptation models discussed in Sections 5.2-5.4. The G × E Interaction module was applied to Data file B on genotype by location cell means for joint regression and AMMI analysis, and GL interaction effects for one-covariate factorial regression. Multi-covariate factorial regression required the performance of multiple regression analyses on the GL effects of the individual genotypes, applying the Regression module with the Multiple data selection option to Data file B. Further calculation for estimation of DF, SS and MS values for insertion in Table 8.2, testing of PC axes by the F_R criterion, model comparison by different procedures etc. was done for original data analysis (as exemplified in Section 8.2 for analysis of transformed data). Winter mean temperature and rainfall amount both reached statistical significance as a covariate in the factorial regression analysis, whereas the selected AMMI model included one PC axis (Table 8.2).

The models ranked in the following order for the ratio of % GL interaction SS to % GL interaction DF: AMMI-1 (ratio = 6.03), joint regression (ratio = 5.80) and two-covariate regression (ratio = 3.98). The criterion based on the sum of the estimated variance of significant GL interaction components of variation ranked the models as follows: two-covariate factorial regression, AMMI-1 and joint regression, with 0.038, 0.031 and 0.030 (t/ha)², respectively (Table 8.2). These indications were fairly inconsistent but suggested placing less emphasis on results of joint regression modelling.

Scores on GL interaction PC axes for genotypes and locations, and expected yield for any genotype-location combination, were outputted by the G × E Interaction module for the selected AMMI-1 model. The graphical expression of genotype adaptive responses in Figure 8.2 (A) (relative to a subset of well-performing cultivars) only required for each genotype the estimation of AMMI-1 nominal yields for the two sites with extreme PC 1 score, i.e. locations 6 and 9 (Table 8.1), and the connection of the two values by a straight line. Nominal yields could be obtained by subtracting the site main effect from the outputted AMMI-1 expected yields, i.e. adding 0.62 t/ha to values for site 6 (the site mean yield and the grand mean being 1.24 and 1.86 t/ha, respectively), and subtracting 1.59 t/ha from values for site 9 (the site mean yield being 3.45 t/ha). Alternatively, these values could be estimated according to formula [5.4], i.e. from values of site PC 1 score (Table 8.1) and genotype mean yield and PC 1 score (Table 8.4), using worksheets with IRRISTAT or other software. For example, the nominal yield of ‘Simeto’ on site 14 is:

2.084 - (0.432 × 0.277) = 1.964 t/ha

AMMI modelling of adaptive responses revealed the occurrence of crossover interaction between top-ranking entries (Fig. 8.2 [A]), which, however, would imply the definition of only two mega-environments (assuming no yield instability of best-yielding material, and two recommended cultivars in all cases):

• one including the high-rainfall, relatively warm sites 1 and 9 (Table 8.1), in which ‘Simeto’ and ‘GTA Dur’ are the best-ranking cultivars;

• the subregion including almost all remaining locations, in which ‘GTA Dur’ and ‘Chen’s’ could be recommended.

Several cultivars were equally well-performing on site 6 (Fig. 8.2 [A]).

FIGURE 8.2. - Nominal yield (A) and nominal yield reliability as lowest yield expected in 80% of cases (B), of eight cultivars modelled as a function of the scores on the first GL interaction PC axis of 14 locations

Source: Annicchiarico et al., 2002a.

The mean value of Dunnett’s one-tailed critical difference (P < 0.20) for entry comparison at a specific site PC 1 score value, calculated as described in Section 5.3 and equal to 0.68 t/ha (from M_err = 0.433 and associated DF = 322, t’ ≈ 2.06, and N = 8), may lead to a wider set of recommended genotypes. Owing to reasons discussed at the end of Section 5.6 and the observed correlation between mean yield and PC 1 score of sites for original yields (r = 0.79, P < 0.01, applying the Regression module on Data file D), the critical difference tends, in this case, to be too large for comparisons at low PC 1 values and too small for comparisons at high PC 1 values of sites. The site score on PC 1 was also correlated positively with winter mean temperature and rainfall, and negatively with altitude.

Genotype regression coefficients and intercept values, and expected GL interaction effects, were outputted by IRRISTAT’s multiple regression analysis of GL effects performed for each genotype for the selected two-covariate factorial regression model. Nominal yield responses cannot be expressed in a graph for this bidimensional model. However, the expected pair of winning genotypes was calculated for each combination of 10 rainfall by 10 winter temperature levels within the range of site values for the covariates (Fig. 8.3). For each combination (represented as a cell in the graph), the nominal yield of each of 10 well-performing genotypes was estimated according to equation [5.7], i.e. from the estimated adaptation parameters of genotypes (Table 8.4) and the relevant pair of covariate values, using worksheets on Data file G (for each genotype, a new variable was created for calculated values). Cells with the same winning genotypes were grouped together. Figure 8.3 reveals the occurrence of crossover interactions between top-ranking entries in the range of site covariates recorded in test years (Table 8.1). It also shows the general interest of ‘GTA Dur’ and the specific interest of ‘Ofanto’, ‘Chen’s’ and ‘Simeto’, depending mainly on the level of winter cold on sites. Landraces and old cultivars (such as ‘Bidi 17’, ‘Polonicum’ and ‘M. Ben Bachir’) and the variety ‘Om Rabi 9’ are preferable only for extremely harsh conditions (Fig. 8.3). There are actually only very few dry, warm sites where ‘Sahel 77’ could be recommended (Table 8.1).

FIGURE 8.3 - Expected pair of top-ranking cultivars for yield at each combination of 10 annual rainfall by 10 mean winter temperature values of locations, according to a two-covariate factorial regression model

Source: Annicchiarico et al., 2002a.

Nominal yields from factorial regression modelling could also be estimated using worksheets for each genotype-location combination as a function of the covariate values recorded in test years, either summing up the genotype mean yield and the expected GL effect stored in Data file B (the GL effects can be stored in the data file as fitted values in the final multiple regression analysis of the selected model performed by the Regression module), or calculating them using equation [5.7]. For example, the estimated nominal yield of ‘Simeto’ on site 14 is:

2.084 - 0.941 + (0.00028 × 297) + (0.1008 × 6.4) = 1.871 t/ha

following the information on relevant parameters reported in Tables 8.1 and 8.4. Calculation confirmed the occurrence of crossover GL interaction between top-ranking entries according to this model.

It is noteworthy that the removal of noise effects by AMMI or factorial regression modelling allowed for a notable simplification of possible recommendation domains. The recommendation based on observed entry values would imply a specific pair of best-ranking genotypes for almost every site (Annicchiarico et al., 2002a).

Finlay and Wilkinson’s regression coefficients of genotypes (outputted from the G × E Interaction module) are also reported in Table 8.4. Going from joint regression to factorial regression through AMMI modelling, the information on specific adaptive responses of genotypes becomes more definite. For example, ‘GTA Dur’ and ‘Simeto’ show: similar responses with respect to site mean yield; slightly different responses according to the genotype PC 1 score; and distinctly different responses with respect to relevant environmental variables, revealing the preference of ‘GTA Dur’ for high-rainfall areas and ‘Simeto’ for mild-winter zones (Table 8.4).

Variation in genotype yield stability

The occurrence of significant GLY interaction in the combined ANOVA (Table 8.2) implies the assessment of, and the comparison among, genotype values of temporal yield stability, to decide whether this characteristic should be taken into account in variety recommendation for each subregion (option 3 vs. 4 in Fig. 2.4). The variance of the average within-site temporal variation of yield values (S_y(l)²) was estimated for each genotype according to formula [7.3], following a separate ANOVA for each genotype performed on plot values and including the factors, location and year within locations (applying the ANOVA module with the Multiple data selection option on Data file A). The S_y(l)² values ranged between 0.510 and 0.167 (t/ha)2 for the 24 genotypes. The ratio of these values (equal to 3.05) was distinctly lower than the critical value of Hartley’s test (greater than 6 for DF = 14, P < 0.05 and 24 compared variances in Table 7.1), suggesting the absence of sizeable variation among genotypes in yield stability. The top-yielding genotype, ‘GTA Dur’, tended to lower Type 4 stability (Table 8.4), but the difference became statistically significant (P < 0.05) according to Fisher’s bilateral test only with respect to the most stable genotype, which was, however, one of the lowest yielding. ‘Om Rabi 9’, the most stable-yielding of the well-performing cultivars (Table 8.4), did not differ from ‘GTA Dur’ even at P < 0.10 (observed variance ratio = 2.02; critical F value for the bilateral test at P < 0.10, coincident with tabular F at P < 0.05, for 14 DF at the numerator and denominator = 2.48). Thus, while genotype differences in temporal yield stability could not be completely ruled out, they appeared modest among well-performing material, suggesting the recommendation for subregions based on the mean value of adaptive responses (option 3 in Fig. 2.4).

Scaling-up of results and final cultivar recommendations for subregions

Recommendations for each of the three sites with only one year’s data could derive from the indications obtained for the most similar location included in the analysis of adaptation, assessing the degree of similarity by the genetic correlation coefficient (see Section 5.7). However, the availability of GIS data was exploited for scaling up spatially (also for these sites) and temporally (for all test sites) the genotype adaptive responses as modelled from data of 14 sites and two years.

For the AMMI-1 model, the first step implied the stepwise multiple regression of the site PC 1 score as a function of mean values of the environmental variables recorded in test years, specifying the forward selection method in the Regression module (applied to Data file D). The dependent variable could accurately be described (R² = 0.81) by the following model, including two variables, selected for a stepwise test level P < 0.10:

PC 1 = -2.025 + (0.00176 × rainfall) + (0.1612 × winter mean temp.)

The second step required the calculation of site PC 1 scores for new sites or test sites based on long-term GIS values of the two environmental variables. For example, the expected PC 1 score for site 14 (rainfall amount = 460 mm, and winter mean temperature = 8.7°C in the long term) is, for an average year, equal to 0.187 (distinctly different from the PC 1 score = -0.277 estimated for test years, which were unusually cold and dry for this site - see Table 8.1). Finally, nominal yields of well-performing genotypes were estimated for all sites by formula [5.4] using the long-term-based estimate of site PC 1 score, and locations having the same pair of best-ranking cultivars were included in the same mega-environment. For example, the nominal yield of ‘Simeto’ on site 14 is now:

2.084 + (0.432 × 0.187) = 2.165 t/ha

As a simpler alternative, the top-ranking material could be identified visually in Figure 8.2 (A) as a function of the predicted PC 1 score of the site. The up-scaling procedure indicated two major subregions based on pairs of top-ranking cultivars, one for ‘GTA Dur’ and ‘Chen’s’ and the other for ‘GTA Dur’ and ‘Simeto’, together with a few additional subregions of very limited extension placed in areas with severe drought and cold stress (data not reported).

The scaling-up of results was more straightforward through the factorial regression modelling, which implied the estimation of nominal yields for well-performing cultivars according to formula [5.7] based on the estimated genotype parameters in Table 8.4 and the long-term values of the two covariates for new sites or test sites. For example, the estimated nominal yield of ‘Simeto’ on site 14 is now:

2.084 - 0.941 + (0.00028 × 460) + (0.1008 × 8.7) = 2.149 t/ha

The results confirmed the presence of three major subregions having ‘GTA Dur’ as a common variety and ‘Chen’s’, ‘Simeto’ and ‘Ofanto’ as specifically recommended cultivars (data not reported). ‘Simeto’ showed a slight increase in the number of test sites where it could be recommended following the temporal scaling-up of responses, because winter temperatures tended towards higher values in the long term than in test years, and this genotype is more adapted than ‘Chen’s’ and ‘Ofanto’ to milder winters (Table 8.4).

The GIS facility allowed for a fine-tuned definition of each subregion as the area on a geographical map including sites with the same pair of expected top-ranking cultivars. Graphs, such as Figure 8.3 for factorial regression or Figure 5.10 for AMMI modelling (also issued from worksheet calculation of expected nominal yields of genotypes as a function of levels of the two environmental variables using Data file G), offer an alternative means for defining recommendations - which can be read as a function of long-term values of climatic variables on each site. For practical use, a denser grid of points for the variables is recommended. AMMI analysis-based indications exploit the above equation for predicting the site PC 1 score, applying it to each combination of levels of climatic variables (see Section 5.8). For combinations in Figure 5.10, the predicted PC 1 score values ranged between -0.94 (for rainfall = 250 mm and temperature = 4°C) and 1.30 (for rainfall = 700 mm and temperature = 13°C). The wider range of negative values compared with that for test year data (Table 8.1) justifies the appearance of some additional recommendation domains (of limited extension, and relative to ‘Om Rabi 9’, ‘M. Ben Bachir’ and ‘Polonicum’ - Fig. 5.10). The inconsistency in recommended genotypes between indications from AMMI (Fig. 5.10) and factorial regression (Fig. 8.3) is larger on dry warm sites on the one hand and wet cold sites on the other (represented on the right upper corner and the left lower corner, respectively, in the figures). Modelling for these conditions, which are scarcely represented on test sites (Table 8.1), was less accurate.

For practical use, AMMI- and factorial regression-based indications that differ for a given area may be integrated for farmer recommendation (contemplating more than two recommended varieties). Data from a third cropping season, used for comparing the different recommendation procedures, have shown a slight advantage of the factorial regression approach over the AMMI-based one, and the distinct advantage of both modelling approaches over the conventional one of recommending, for each test location, the best-ranking cultivars on the basis of observed yield data (Annicchiarico et al., unpublished results).

Exemplification of cultivar recommendations over the region

Procedures for cultivar recommendations that are relative to options 1, 2 and 4 in Figure 2.4 (not useful in this case) are briefly exemplified with respect to the same data set. If modelling of adaptation patterns had shown no crossover GL interaction of top-ranking genotypes, the possible choice between options 1 and 2 in the presence of significant GLY interaction would have depended on the occurrence of yield instability for some high-yielding entry (Fig. 2.4).

Yield stability as environmental variance (S²) and mean yield of genotypes could conveniently be estimated from data of all available test environments (i.e. 31) for recommendation over the region. Environmental variance values ranged between 1.23 and 0.44 (t/ha)² for the 24 genotypes. The least stable-yielding entry (‘GTA Dur’) was also the top-ranking for mean yield (Table 8.4). Of the relatively well-performing cultivars, the most stable was ‘Om Rabi 9’ (S² = 0.76 [t/ha]²), which was significantly more stable than ‘GTA Dur’ at P < 0.01, according to Ekbohm’s test. In the test, the observed variance ratio = 1.23/0.76 = 1.62 was compared with the critical F value ≈ 1.36, relative to tabular F at P < 0.005 for the bilateral test and DF = 307 for the numerator and denominator. The DF value was calculated as described in Section 7.1, after estimating the correlation r = 0.953 between yield values of the two genotypes across the environments (through the Regression module applied to a novel data file reporting the matrix of genotype-environment cell means with environments as rows and genotypes as columns). This result supports cultivar recommendation over the region based on Kataoka’s measure of yield reliability. Other well-performing genotypes were significantly (P < 0.05) more stable-yielding than ‘GTA Dur’ (Table 8.4).

It should be noted that Type 1 and Type 4 stability of the 24 genotypes (expressed as square root of S² and S_y(l)² values, respectively, to approximate their distribution to the normal one) were not correlated (r = 0.24), highlighting the bearing of the inclusion or exclusion of responses to locations on the results. There was no definite relationship between temporal stability and consistency of response to locations, as shown by the lack of correlation (r = 0.26) between square root of S_y(l)² values and absolute value of GL interaction PC 1 score (score near zero indicating consistent response).

In this example, Kataoka’s index of reliability estimated the lowest yield that can be expected in 80 percent of cases (i.e. P = 0.80). The index value was calculated from the mean yield across environments and the square root of the S² value of each genotype according to formula [7.4]. Mean values of genotypes are reported in Table 8.4 (where they are compared through the critical difference d = 0.185 t/ha for Dunnett’s one-tailed test at P < 0.20, calculated after performing a combined ANOVA according to Model 2 in Table 4.1 using the ANOVA module on Data file C). For example, the value I = 1.250 for ‘GTA Dur’ derives from the equation: 2.182 - (0.840 × 1.109), Z_(P) being 0.840 for P = 0.80. The stable-yielding genotypes ‘Om Rabi 9’ and ‘Ofanto’, which tended towards lower merit than ‘GTA Dur’ and ‘Chen’s’ on a mean yield basis, reached a similar merit in terms of yield reliability (Table 8.4). Dunnett’s critical difference for separation of genotype mean values across environments (d = 0.185 t/ha) may also provide a rough indication of significantly less reliable entries compared with the top-ranking one (i.e. those where I < 1.10 in Table 8.4). In any case, a set of recommended cultivars is likely to include no more than the four or five best-yielding genotypes.

For extension to farmers, the current estimates of reliability (relating to yield values averaged across environments) may be expressed in terms of percentage values or yield differences relative to the most reliable cultivar. For example, the lowest yield provided by ‘Simeto’ in four cases out of five is about 10 percent lower than that provided by ‘Chen’s’ (Table 8.4).

Exemplification of cultivar recommendations for subregions based on yield reliability

A final possibility is represented by specific recommendations for subregions based on the reliability of nominal yield responses (option 4 in Fig. 2.4). The lowest nominal yield expected in 80 percent of cases is reported in Figure 8.2 (B) for AMMI-modelled adaptive responses of a set of well-performing genotypes. The only difference in the calculation compared to mean values of nominal yield responses (Fig. 8.2 [A]) is represented by the addition of the term: -(Z_(P) S_y(l)), where Z_(P) = 0.84, and S_y(l) = square root of the genotype S_y(l)² value reported in Table 8.4. For example, the reliability of nominal yield of ‘Simeto’ on site 14, based on the PC 1 score of the site in test years (Table 8.1) and the estimated genotype parameters (Table 8.4), is:

2.084 - (0.432 × 0.277) - (0.840 × 0.648) = 1.420 t/ha

For each genotype, the response in Figure 8.2 (B) may be represented by connecting through a straight line the values calculated for the two sites with extreme PC 1 scores. Compared with results for nominal yield, those for nominal yield reliability show a relatively better response of cultivars tending towards higher temporal yield stability, such as ‘Om Rabi 9’ (top-ranking in the stressful site 6) and ‘Bidi/Waha/Bidi’, and a relatively worse response of the somewhat less stable cultivar ‘GTA Dur’ (Fig. 8.2).

Likewise, the estimation of genotype responses for reliability of nominal yield according to a factorial regression model requires the addition of the term -(Z_(P) S_y(l)) to previous formulae. For example, the lowest nominal yield of ‘Simeto’ expected in 80 percent of cases on site 14 may be calculated as follows on the basis of test year values of the two covariates in the model (Table 8.1) and the estimated genotype parameters (Table 8.4):

2.084 - 0.941 + (0.00028 × 297) + (0.1008 × 6.4) - (0.840 × 0.648) = 1.327 t/ha

For final recommendations, results can be scaled up as a function of long-term values of the two climatic variables on new sites or test sites. As in earlier up-scaling procedures, the information on climatic variables can be inputted directly in the equations for factorial regression modelling and indirectly in those for AMMI modelling (calculating, in the latter case, the expected site PC 1 score as a function of climatic variables). The only difference compared to previous formulae is the inclusion of the term (Z_(P) S_y(l)), in order to estimate the reliability of nominal yields. Procedures that imply the AMMI- or factorial regression-based estimation of reliability values of genotypes for each individual site can conveniently be resumed by a geographical map reporting subregions (as areas with the same pair of most-reliable genotypes). Figure. 8.2 (B) could be used in the same context by reading the most reliable cultivars as a function of the expected site PC 1 score.

Other up-scaling procedures contemplate the AMMI- or factorial regression-based estimation of nominal yield reliability of genotypes for each combination of levels of rainfall and winter mean temperature, reporting the most reliable genotypes for each combination in a graph (in which cultivar recommendations can be read as a function of long-term values of the climatic variables). Calculations can be limited to a set of well-performing genotypes using worksheets on Data file G.

For example, the AMMI-modelled nominal yield reliability of ‘Simeto’ at rainfall = 450 mm and winter mean temperature = 9°C (based on the above general equation for prediction of site PC 1 score and the estimated parameters of the genotype in Table 8.4) is:

2.084 + 0.432 [-2.025 + (0.00176 × 450) + (0.1612 × 9)] - (0.840 × 0.648) = 1.634 t/ha

while the factorial regression-modelled value for the same genotype and conditions is:

2.084 - 0.941 + (0.00028 × 450) + (0.1008 × 9) - (0.840 × 0.648) = 1.632 t/ha

The most reliable pair of genotypes for each combination of climatic variables according to the factorial regression model is reported in Figure 8.4. The comparison with indications for mean adaptive responses based on the same model (Fig. 8.3) confirms the larger extent of conditions under which ‘Om Rabi 9’ and ‘Bidi/Waha/Bidi’ could be preferable, were temporal yield stability taken into account.

Values of nominal yield reliability outputted for a set of cultivars as a function of inputted long-term values of climatic variables (as may be the case using a simple Decision-aid System) may be expressed in terms of yield differences relative to the most reliable genotype, to facilitate the extension of results.

FIGURE 8.4 - Expected pair of top-ranking cultivars for yield reliability at each combination of 10 annual rainfall by 10 mean winter temperature values of locations, according to a two-covariate factorial regression model (yield reliability as lowest yield expected in 80% of cases)