Introduction

Sugarcane (Saccharum spp.) is the world’s most important biomass-producing crop. It is a cash crop grown over semiperennial durations, with harvesting occurring at yearly intervals from time of planting. Breeding programs for the cane industry aim to select new genotypes with local or broad adaptation to provide high yields over several crop years within target environments representative of different cultivation areas. Breeding programs consist of several successive stages of selection using clonal multiplication. In the most advanced selection stages, the best elite candidates are tested in multienvironment trials (MET) at representative locations using replicated experiments. MET carried out over several crop years allow testing of genotype × location (GL) and genotype × crop year (GC) interactions. These two components of the genotype × environment (GE) interaction determine the yield potential and yield stability of sugarcane varieties. They have been quantified for cane yield and sucrose content or for derived sucrose yield in breeding programs at various locations, including Australia (Jackson and Hogarth 1992; Mirzawan et al. 1994), Florida (Milligan et al. 1996; Glaz and Kang 2008), South Africa (Parfitt 2000; Ramburan et al. 2012a), Venezuela (Rea and De Souza Vieira 2002), and Argentina (Mariotti and Clariana 1994). All these examples of MET studies encompassed locations spread over relatively large cultivated areas, necessarily made up of environments that are more or less heterogeneous in various regards, such as soil type (Milligan et al. 1996; Glaz and Kang 2008) or cultivation practices (Mirzawan et al. 1994). Studies providing formal statistical tests of interactions revealed that both GL and GC effects were always significant for all traits surveyed, except on one occasion (Parfitt 2000). Moreover, GL interactions were more important than GC interactions in these reports.

G × E interactions may cause changes in the relative ranking of genotypes across sites and/or crop years in trials and complicate identification of superior cultivars by confounding determination of true genetic values. When G × E interactions exist, their statistical significance and precise characteristics must be investigated in detail to assess implications for selection strategies and help optimize resource allocation across locations and years. As a first step, analysis of variance of MET data provides a general picture of the influence of the different factors underlying phenotypic variation. In a second step, genotype plus genotype × environment interaction (GGE) biplot analysis (Yan and Tinker 2005, 2006) is a very useful graphical tool to investigate in detail the relationships between environments and the pattern of the response of genotypes across environments. This popular visualization technique for MET data has been used in several sugarcane programs to (i) investigate the similarity of environments and their ability to discriminate genotypes (Glaz and Kang 2008; Ramburan et al. 2012a; Luo et al. 2015), (ii) identify redundant sites or megaenvironments and analyze the stability of the site response across series of genotypes (Ramburan et al. 2012a, b), and (iii) visualize the performance rank and stability of genotypes across environments for the purposes of decision-making regarding release of new cultivars (Glaz and Kang 2008; Shandu et al. 2012; Klomsa-ard et al. 2013; Luo et al. 2015).

On Réunion Island, production of sugarcane is scattered over many different agroclimatic regions on both leeward (dry) and windward (wet) coasts. These different regions range in altitude from sea level to highlands up to 900 m and consist of vastly different soil types in terms of physical properties and chemical fertility. Commercial sugarcane genotypes are selected and released by eRcane research institute. eRcane currently operates seven selection programs (eRcane 2009), strategically located in the major sugarcane growing areas (Table 1). Programs run concurrently to identify promising genotypes suited for each zone and currently take 14 years from time of initial cross to new cultivar release. During the last four testing years, these seven programs have shared a common final selection stage, gathering a set of elite genotypes selected from each of the seven selection sites. Each year, a new MET series of genotypes is planted across the seven locations and tested over several crop years. Since the inception of the present-day MET network in the late 2000s (eRcane 2009), no study has been conducted to evaluate the GE interactions and the relationships among the locations in terms of genotypic response.

Table 1 Climate and soil description of the seven eRcane locations used to test advanced MET selection series

The present study aims to investigate the characteristics of the present-day MET selection program on Réunion Island by retrospective analysis of some recent genotype series tested across the whole network of selection stations. The objectives are to: (i) assess variance components of the major quantitative traits relative to yield components and estimate the importance of GE interactions, (ii) study relationships existing between environments in terms of genotypic response, (iii) test the added value that the GGE biplot statistical tool can bring to provide support in decision-making for selection, and (iv) identify potential areas for optimizing the current MET program.

Materials and methods

Multienvironment trial (MET) dataset

The MET analyzed in this study were carried out at a network of seven selection sites on Réunion Island, France: La Mare (LM), Saint-Benoit (SB), Menciol (MN), Étang-Salé (ES), Le Gol (GL), Vue-Belle (VB), and Saint-Philippe (SP). These sites cover a wide range of ecologies of production representative of the main sugarcane growing areas. The main climatic and soil type characteristics of these seven locations are summarized in Table 1. Sites represent either coastal zones (LM, SB, ES, and GL) or zones at low (SP) or medium (MN) altitude, as well as highland (VB) characterized by cooler temperatures. The variability in terms of soil type and physical characteristics for cane cultivation is relatively large. Chemical fertility is also highly variable, being either favorable (ES, GL), satisfactory (LM, SB, SP), or very low (VB, MN).

The MET dataset consists of four consecutive genotype series (S00, S03, S04, and S05) planted in the years 2011–2014 at the seven sites. Each MET series was tested in a randomized complete block design with four replications at all locations. Each plot had area of 45 m2 with three rows, with length of 10 m and width of 4.5 m. As indicated in Table 2, MET series were phenotyped and harvested at all sites for between one and three crop years, according to the year of their planting. Due to constraints in terms of land resources and the high cost of the whole breeding program, the material composing a final MET series usually represents a set of genotypes that is not balanced across all sites. On average, only about half of the clones tested at one location were also tested at the other six locations. Altogether, the four MET series studied herein represent a balanced number of 47 genotypes tested across all sites. Depending on the site considered, two (LM, MN, SP, VB, SB), three (ES), or four (GL) of the best local cultivars were also planted as common standards in the trials of the four genotype series (see footnote b of Table 2).

Table 2 Summary of the four advanced MET selection series (S00, S03, S04, and S05) tested at the seven locations described in Table 1

At the end of each crop year, four quantitative traits were recorded at individual plot level: tonnes of cane per hectare (TCH, Mg ha−1), fiber content as percentage of fresh weight (FIB, %) estimable recoverable sugar (ERS, %) and an economic index (EI). All millable stalks from each variety plot were manually cut and weighed using a digital scale mounted on a tractor-operated hydraulic boom. TCH was determined from plot weights divided by plot areas. From each plot, a sample of 18 randomly selected stalks was used to determine ERS and FIB at eRcane laboratory using the standard hydraulic press method (Hoarau 1969): FIB was determined from the weight of the press cake from a 500-g shredded subsample. ERS was calculated from the FIB value and Brix (Bellingham RFM340 refractometer) and Pol (Polaser SR 64 polarimeter) measurements of the extracted juice according to conventional calculations used in the local cane industry adapted from Saranin (1986). The conventional economic index (EI) used by eRcane to rank the cultivation merit of candidate genotypes was calculated from the TCH and ERS using the formula EI = TCH × (ERS − 4), being indicative of farming profit net of production costs (Hugot 1958).

Variance components analysis

Variance components for each trait were assessed using mixed linear models designed to estimate GL and GC interaction effects. This was first done separately for each MET series consisting of unbalanced numbers of genotypes among sites, with the following mixed linear model fit to each trait:

$$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{Y}_{jklm} = \mu + L_{j} + C_{k(j)} + R_{l(j)} + \underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{G}_{m} + \underline{GL}_{mj} + \underline{GC}_{mk} + \underline{GLC}_{mjk} + \underline{e}_{jklm} \quad ({\text{model}}\;1),$$

where Y jklm is the observation of genotype m in replicate l in crop year k at location j, μ is the grand mean, L j is the location main effect at j, C k(j) is the crop year effect for k at location j, \(R_{l(j)}\) is the replication effect for l at location j, G m is the main effect for genotype m, GL mj is the effect of the genotype × location interaction, GC mk is the effect of the genotype × crop year interaction, GLC mjk is the effect of the genotype × crop year × location interaction, and e jklm is the error. Effects considered as random in the model are underlined, while others were fixed. In the particular cases of series S04 and S05, which were studied in a single year, the crop year effect (C k(j)) was dropped from the model.

Combined global analysis of all four MET series was carried out using the subset of 47 genotypes that were commonly tested at all seven sites using the following mixed linear model fit to each trait:

$$\underline{Y}_{ijklm} = \mu + S_{i} + L_{j} + C_{k(i,j)} + R_{l(i,j)} + \underline{G}_{m} + \underline{GL}_{mj} + \underline{GC}_{mk} + \underline{GLC}_{mjk} + \underline{e}_{ijklm} \quad ({\text{model}}\,2),$$

where S i is the main effect of series i, while the other terms and indices are the same as in model 1.

All models were performed using the MIXED procedure of SAS 9.3 (SAS Institute, Cary, NC), and variance component estimates were evaluated using the restricted maximum likelihood (REML) procedure. Use of the COVTEST option in the model statement provided standard errors of variance components and a derived Wald Z-test of their statistical significance (Littell et al. 2006). In both models 1 and 2, measurements acquired on the same individual plot across successive crop years in series S00 and S03 were considered as repeated measurements (longitudinal data).

When analyzing series S00 and S03 with model 1, four R matrices of the variance–covariance (VCV) of random error terms were tested to select the structure that modeled the data optimally: (i) using the UN option to allow all VCV parameters to be unstructured (pairs of within-genotypes errors having their own correlations), (ii) using the UN(1) option to specify the constraint of covariances of error terms between crop years to be null, (iii) using the first-order autoregressive AR(1) option with the constraint of correlation between adjacent within-genotypes errors |ρ| < 1, and (iv) using a component symmetry (CS) option having the constraint of constant variances and constant covariances between error terms. The goodness of fit values for the UN, UN(1), AR(1), and CS models were compared using the Akaike information criterion (AIC) (Akaike 1974), which balances model fit versus number of parameters.

When analyzing all four MET series together, an UN structure was considered for the R matrix.

Broad-sense heritability H or genetic repeatability refers to the extent to which the phenotype is determined by its genotype (Falconer and Mackay 1996). At MET level, H was calculated for each trait with variance components estimated from each model, as a ratio of the genetic variance (\(\sigma_{\text{G}}^{2}\)) to phenotypic variance of genotype means (\(\sigma_{\text{P}}^{2}\)):

$$H = \frac{{\sigma_{\text{G}}^{2} }}{{\sigma_{\text{P}}^{2} }},$$

where \(\sigma_{\text{P}}^{2} = \sigma_{\text{G}}^{2} + \frac{{\sigma_{\text{GL}}^{2} }}{L} + \frac{{\sigma_{\text{GC}}^{2} }}{C} + \frac{{\sigma_{\text{GLC}}^{2} }}{LC} + \frac{{\sigma_{\text{e}}^{2} }}{LCR},\) with \(\sigma_{\text{GL}}^{2}\), \(\sigma_{\text{GC}}^{2}\), and \(\sigma_{\text{GLC}}^{2}\) designating GL, GC, and GLC interaction variances, respectively, and considering L = 7 locations, R = 4 replicates in both models, with C = 1–3 crop years in model 1 depending on the MET series considered and C equal to the harmonic mean of the number of crop years of the different series in model 2 (Holland et al. 2003, pp. 64–65).

Model 2, which was used to analyze the combined MET series, was considered to provide more representative estimates of variance components than those obtained from model 1, because of the higher number of genotypes used to compute the estimates (47 instead of about 20) and the subsequent smaller risk of high individual influence of genotypes on variance component partitioning. Using the variance component estimates obtained from model 2, we simulated the broad-sense heritability (H) of traits at MET level using various combinations of number of locations (L) and crop years (C) up to the maximum numbers in the experiments (L = 7 and C = 3) and considering a standard trial design with R = 4 replications. The changes in H when the number of locations and replications were increased were plotted graphically to explore the effect of these changes and assess the influence of experimental effort on the accuracy of the mean performance of genotypes across all sites.

GGE biplots

To visualize the relationships existing among selection sites and the performance of some of the most promising candidate genotypes, GGE biplots were produced for the 21 genotypes of both series S00 and S03 (10 and 11 genotypes, respectively) that were tested across all sites (Table 2). To this end, we restricted the trait data analyzed to the first two crop years in order to consider space–time information perfectly balanced (7 locations × 4 replications × 2 crop years × 21 genotypes) to visualize relationships among environments and genotypes without risks of bias. In a first step, using cultivar checks common to both series within each site (see footnote b of Table 2), both series were analyzed at each site using the following linear mixed model:

$$\underline{Y}_{iklm} = \mu + \underline{S}_{i} + \underline{C}_{k(i)} + \underline{R}_{l(i)} + G_{m} + \underline{e}_{iklm} \quad ({\text{model}}\,3),$$

where the symbols and indices designate the same effects as in model 2 and in which the genotype effect was fixed and other effects were random. At each location j, the adjusted mean (\(\bar{Y}_{mj}\)) of genotype m (free of any trial/series membership effects) was inferred by adding to the best linear unbiased estimate (BLUE) of genotype the overall mean (µ) of the model (Littell et al. 2006, p. 211). Pearson’s correlation coefficients between traits and across sites were calculated using average adjusted means of genotypes for the seven sites. For each trait, the contribution of each genotype m to the genotype × location interaction was estimated using Wricke’s (1965) ecovalence stability index (W m ), and the contribution of each location j to the interaction was estimated using a similar index (W j ) by interchanging locations and genotypes as follows:

$$W_{m} = \mathop \sum \limits_{j = 1}^{7} (\bar{Y}_{mj} - \bar{Y}_{m .} - \bar{Y}_{.j} - \bar{Y}_{..} )^{2},$$
$$W_{j} = \mathop \sum \limits_{m = 1}^{21} (\bar{Y}_{mj} - \bar{Y}_{m .} - \bar{Y}_{.j} - \bar{Y}_{..} )^{2},$$

where \(\bar{Y}_{m .}\), \(\bar{Y}_{.j}\), and \(\bar{Y}_{..}\) are the adjusted means of genotypes and locations, and overall mean, respectively. The two-way data table of adjusted genotype means × locations (\(\bar{Y}_{mj}\)) was then centered to the mean trait value (µ j ) of each environment j and divided by its standard deviation (s j ) to obtained a “standardized GGE matrix” of the genotype main effect (G) and genotype × environment interaction (GE). This standardized GGE matrix was subjected to singular value (SV) partitioning between the genotype and environment eigenvectors using the general models of Yan (2002) and Yan and Tinker (2005, 2006):

$$(\bar{Y}_{mj} - \mu_{j} )/s_{j} = \mathop \sum \limits_{p = 1}^{n} \lambda_{p}^{{f_{p} }} \alpha_{mp} \lambda_{p}^{{1 - f_{p} }} \gamma_{jp} + e_{mj},$$

where λ p is the SV of the pth principal component (PC), α mp and γ jp are the respective eigenvectors of genotype m and environment j for PC p, e mj is the residual associated with genotype m in environment j, and f p is a partition factor for PC p, equal to 0.5 for all PCs p, to visualize genotype scores and environmental scores in the same units for both PC1 and PC2 (Yan 2002).

GGE biplots provide, at a glance, the ranking of all genotypes regarding their performance in any environment, when visualizing the distribution of their positions when projected orthogonally onto each environment axis (Yan and Tinker 2006). To assess the efficiency of our GGE biplots for visualizing the trait performance of genotypes in each environment, we calculated for each environment a coefficient of correlation between the rank of genotypes in the original data and their apparent rank along the environment axis. The latter was precisely inferred by computing genotype abscissa on the axis of each environment j resulting from the scaler product between the vector of each genotype m (\(\vec{u}_{m}\)) and the environment vector (\(\vec{v}_{m}\)) divided by the length of the latter (\(\vec{u}_{m} \cdot \vec{v}_{m} / \, \| \vec{v}_{m} \|\)) (Spiegel et al. 2009).

Results

Variance components analysis

The AIC values for model 1 associated with the different R structures tested in trait analyses for series S00 and S03 are presented in Table 3. In both series, the TCH and EI data were best modeled using the UN structure and the FIB data when using the UN(1) structure (lowest AIC values), while the ERS data were best modeled when using either the UN (S00) or UN(1) (S03) structure.

Table 3 Number of fit parameters and Akaike information criteria (AIC) of different R matrices for model 1 applied to S00 and S03 MET data for tonnes of cane per hectare (TCH), estimable recoverable sugar (ERS), fiber content (FIB), and economic index (EI)

Accordingly, Table 4 presents the estimates of variance components obtained from the separate analyses of the four series (model 1) or from the “global” analysis performed on a combination of balanced subsets of genotypes from each series across all sites (model 2). As a general rule, the statistical significance of genotype (\(\sigma_{\text{G}}^{2}\)) and genotype × location (\(\sigma_{\text{GL}}^{2}\)) variance components may fluctuate for some traits across the four individual series (S00, S03, S04, and S05), as well as the genotype × crop year (\(\sigma_{\text{GC}}^{2}\)) and genotype × location × crop year (\(\sigma_{\text{GLC}}^{2}\)) variance components across the two oldest series (S00 and S03). However, estimates of variance components obtained from all these individual analyses (model 1) were convergent with estimates obtained from the “global” model (model 2). In this reference second model, the \(\sigma_{\text{G}}^{2}\) and \(\sigma_{\text{GL}}^{2}\) variance components were always highly significant (P < 0.01 or P < 0.001) for all traits, as well as the \(\sigma_{\text{GLC}}^{2}\) component, except for FIB since this component was estimated to be null. For \(\sigma_{\text{GC}}^{2}\), this component was more or less significant for TCH (P < 0.01) as well as for FIB and EI (P < 0.05), but was not significant for ERS (P > 0.05). Apart from residual variance, the order of importance of the variance components in terms of percentage of phenotypic variance explained was G > GL > GLC > GC for TCH and ERS, G > GL > GC > GLC for FIB, and GL > G>GLC > GC for EI. The genetic variance (\(\sigma_{\text{G}}^{2}\)) represented 31, 24, 23, and only 16%, for FIB, ERS, TCH, and EI, respectively.

Table 4 Estimates of variance components, percentage of phenotypic variance explained, and broad-sense heritability at experimental level for tonnes of cane per hectare (TCH), estimable recoverable sugar (ERS), fiber content (FIB), and economic index (EI) from individual (model 1) or partially global (model 2) analyses of MET series data

Broad-sense heritability with different locations × crop years entry-mean bases

Simulated broad-sense heritability (H) at the maximal experimental level of seven locations and three years reached 0.76 for EI, 0.85 for TCH, 0.88 for ERS, and 0.92 for FIB (Fig. 1). Accordingly, lower entry-mean heritability (H) simulated with fewer locations and years (L ≤ 7 and C ≤ 3) always showed the same ranking of H values among traits, i.e., EI ≪ TCH ~ ERS ≪ FIB. When considering a single location and single crop year (L = C = 1), H was very low for EI (0.26), modest for TCH and ERS (0.38 or 0.41), and relatively high for FIB (0.55). An increase of one or two unit(s) of either L or C (or both) rapidly enhanced the H values of all traits. For L + C ≥ 3, H evolved more slowly toward final plateaus that conserved the initial ranking between traits regardless of which entry-mean basis was examined. The most predictable traits (FIB, ERS, and TCH) required less L + C to attain reliable H estimates compared with the less predictable trait (EI). Location generated a slightly greater increase of H compared with crop year: depending on the trait considered, an increase of two units in location (L = 3 compared with L = 1) resulted in an increase in the range of 0.20–0.26 in H but only 0.08–0.14 when considering a similar increase in crop years (C = 3 compared with C = 1).

Fig. 1
figure 1

Simulated broad-sense heritability (H) of tonnes of cane per hectare (TCH), estimable recoverable sugar (ERS), fiber content (FIB), and economic index (EI) expressed as a function of their variance components (model 2) and combinations of various numbers of locations and crop years

Statistical analysis of GGE data

The adjusted means of the 21 genotypes of the two oldest series (S00 and S03) which were tested in common at the seven locations across two crop years (Tables 5, 6) represented the balanced GGE data further used to construct biplots representative of genotype × environment (GE) relations. This genotype set, which gathered elite candidates previously selected from one of the seven contrasted locations, logically exhibited (i) a relatively high 20.6% coefficient of variation (CV) mean between genotypes for TCH, (ii) more moderate CV means for ERS (7.8%) and FIB (8.9%), and (iii) a subsequent relatively high CV of 20.9% for the derived EI trait. The adjusted means of the 21 genotypes for all locations were subjected to analyses of variance and the trait means for the locations were compared using a Duncan’s multiple-range test at P = 0.05. ANOVA revealed a highly significant (P < 0.0001) variance for location for each of the four yield-related traits (TCH, ERS, FIB, and EI). Statistical ranking of location means varied a lot between TCH (SB > ES ≥ LM ≥ GL > MN = SP = VB), ERS (LM > GL = SP = VB > MN = SB > ES), FIB (GL > MN ≥ LM ≥ ES = VB > SB > SP), and EI (LM ≥ SB ≥ GL > ES = VB > MN = SP) and did not reveal any broad similarity between any two traits. Correlations of traits across locations revealed in the set of genotypes surveyed are presented in Table 7. There was no correlation between TCH and FIB or between ERS and EI, high positive correlation between TCH and EI (0.746), high negative correlation between TCH and ERS (−0.626), and moderate negative correlations between ERS and FIB (−0.367) and between FIB and EI (−0.342). The ecovalence of the seven locations (\(W_{j}\)) varied around its average value (14% = 100%/7 locations) with larger range (rg) for TCH (rg = 20.9%) and ERS (rg = 24.9%) compared with FIB (rg = 12.9%) or EI (rg = 13.4%). For TCH, the smallest contribution to genotype × location interaction was due to location LM (W LM = 5.3%) and the highest to location SB (W SB = 26.2). For ERS, the smallest contribution to the interaction was due to ES (W ES = 6.4%) and GL (W GL = 6.6%) and the highest to MC (W MC = 31.3%). Similarly, for each trait, the individual contribution of genotypes to the genotype × location interaction was variable around the average ecovalence value of W m  = 4.76% (100%/21 genotypes). However, the individual contribution of genotype remained rather modest, since in 85% (77/84) of all cases (four traits considered together) W m did not exceed 7% and in the remaining 15% of cases the highest genotype ecovalence reached a maximum of 17.7% (genotype G7 for FIB).

Table 5 Adjusted mean for tonnes of cane per hectare (TCH) and estimable recoverable sugar (ERS) over the two first crop years of the 21 genotypes of series S00 and S03 tested at the seven locations along with location or genotype ecovalence indexes and variation between genotypes
Table 6 Adjusted mean for fiber content (FIB) and economic index (EI) over the two first crop years of the 21 genotypes of series S00 and S03 tested at the seven locations along with location or genotype ecovalence indexes and variation between genotypes
Table 7 Global phenotypic correlations between the four yield-related traits computed on the mean value across locations of the balanced genotype dataset used for the GGE study

Overview of the distribution of genotypes in GGE biplots

The symmetric scaling of genotype and environment scores in GGE biplots permits direct visualization of the magnitude of genotypic and environmental variations in the same units for both PC1 and PC2 (Yan 2002). Figure 2 presents GGE biplots of the four traits. As a general rule, for all traits, environment vectors always had positive abscissa and the variation among them was first discriminated by PC2 while variation among genotypes was first discriminated by PC1. These first two principal components (PCs) of biplots explained 76.52, 71.55, 90.23, and 63.41% of the total GGE variation of TCH, ERS, FIB, and EI, respectively (Fig. 2). This suggests that a biplot represented by both PC1 and PC2: (i) adequately approximates the GGE data of TCH and ERS, (ii) represents very accurately the GGE data of FIB, but (iii) represents less efficiently the GGE data of EI. For this latter trait, data variability can be represented with a level roughly similar to the other three biplots when taking into account the additional contribution of its PC3 (15.16%). The polygon formed by connecting the genotypes that are further away from the biplot origin contains all genotypes. The wider spread of genotypes over the whole biplot plan for both TCH and EI as opposed to both ERS and FIB reflects the much higher contribution of the respective PC2 for the two former characters (25.58 and 31.23%) compared with the latter two (15.75 and 4.71%) in the variation of their respective GGE data. Comparing the distribution pattern of genotypes in the biplots of TCH and ERS showed abscissa of opposite sign for a very large number of genotypes (17/21) between the two traits. The high frequency of this opposite direction of genotype abscissa between TCH and ERS is in agreement with (i) the negative correlation coefficient between these two traits (−0.626), and (ii) the similarity of the weight of their PC1 axes (50.94 or 55.80%) in the variation of GGE data. Finally, the coordinates of genotypes in TCH and EI biplots were of identical sign for a majority of 14 of the 21 abscissa (8 negative and 6 positive) and a majority of 16 of the 21 ordinates (8 positive and 8 negative). These double findings are perfectly congruent with (i) the positive correlation observed between these two traits (0.746) and (ii) the similar order of magnitude of the contribution of both their PC1 (50.94 and 32.18%) and PC2 (25.58 and 31.23%) axes in the variation of their GGE data.

Fig. 2
figure 2

GGE biplots of the four yield-related traits (TCH, ERS, FIB, and EI) based on data of 21 genotypes of series S00 and S03 tested at seven locations. PC1 and PC2 are the principal component scores on the first and second axis, respectively. The variation accounted for by the axes is shown in brackets. Biplots are based on symmetric scaling between genotypes and environments (scores with the same units for both PC1 and PC2). Environments are indicated by two letters, and genotypes (G) by numbers. The reference frame indicated by dashed arrows represents an average environment vector and its orthogonal vector

The which-won-where pattern

The format of the polygons surrounding all genotypes displays the which-won-where pattern and hence is a succinct summary of the GE pattern of the MET dataset. The rays that are perpendicular to the sides of the polygon (or their extensions) divide the TCH, ERS, FIB, and EI biplots into six, six, six, and nine sectors, respectively, allowing immediate visualization of the similarity between locations for the genotype response. An interesting feature of this view of a GGE biplot is that the genotype on the vertex for each sector had the highest trait value in all environments that fall in the sector, providing that the percentage of GGE data explained by the biplot is high enough to accurately reflect the original data (Yan et al. 2000; Yan 2002). For instance, for the TCH biplot that explained 76.52% of the GGE data, the seven environments fell in one sector for three of them (SP, GL, and ES). The four remaining environments (MN, LM, VB, and SB) fell in the second sector. The vertex genotypes for these two sectors were G19 and G9, respectively, suggesting that these two genotypes had the highest tonnage or nearly the highest in all environments that fell in each respective sector. In fact, both G19 and G9 were by far the best yielding genotypes in all locations falling in their respective sectors (Table 5), except in VB where G9 was the second best genotype just after G10 (with only 1.1 Mg ha−1 difference). The reliability of the which-won-where pattern revealed in biplots by vertex genotypes depends on (i) the cumulated weights of PC1 and PC2 axes, and (ii) the narrowness of the angles that encompass all the vectors of environments falling in each sector of interest. For example, in the biplot of FIB, the weight of the first two axes was rather high (90.23%) and both angles in the two sectors that grouped either four or three locations were each relatively acute. Congruently, the which-won-where pattern identified the genuine highest variety in five environments (G12 in ES, GL, and SB; G20 in MN and SP) but the second highest one in two environments (G12 in LM; G20 in VB). In the biplot of ERS, the first two axes cumulated a weight (71.55%), slightly smaller than that of the FIB biplot, and the angle that encompassed almost all environment vectors (6/7) in a single sector was slightly wider than the two angles in the FIB biplot. Congruently, the which-won-where pattern identified the genuine highest variety in five environments (G2 in LM, SB, SP, and VB; G14 in MN) but the third best in one environment (G2 in GL), and the seventh best in the remaining environment (G2 in ES). In the biplot of EI, the total weight of both first axes dropped to 63.41% of the total GGE variation and the angle that encompassed four of the seven environment vectors in the same sector (GL, ES, LM, and SP) was significantly wider than the angles examined in all of the other three biplots. Congruently, the which-won-where pattern identified the genuine best variety only in three environments (G19 in LM and SP; G10 in VB) but the second best in one environment (G19 in ES), the third best in one environment (19 in GL), and the fourth best in the two remaining environments (G10 in MN; G9 in SB). To obtain a more global view of the ability of biplots to display reliable classifications of the performance of all genotypes, we calculated a correlation coefficient between the genotype ranking along each environment axis and the corresponding genotype ranking in the initial data (Tables 5, 6) and computed an overall mean across all seven environments. This mean correlation coefficient for genotype ranking reached 0.93 for FIB, 0.84 for both TCH and ERS, and 0.50 for EI. The hierarchy of these four coefficients of correlation (FIB > TCH = ERS > EI) and the magnitude of their relative differences are perfectly in alignment with the hierarchy and differences observed in the percentage of variation, respectively, explained by the biplots of the four traits (90.23% > 76.52% ≈ 71.55% > 63.41%).

Interrelationship among environments

Correlation coefficients among the seven environments are presented in Tables 8 and 9 with bold characters indicating values that are statistically different from zero (P < 0.05). The vector view of GGE biplots provides a succinct summary of the interrelationships among environments. For instance, in the biplot of TCH, the seven environments fell in two sectors that grouped either three or four environments. More or less acute angles between vectors could only be observed between environments within a same sector, congruently with the facts that (i) all correlation coefficients between any two locations within each sector were significant, and (ii) conversely correlations between any two environments belonging to different sectors were not significant, with the sole exception of the correlation between MN and ES [which can easily be explained by the coordinates of these two environments on PC3 with significant value and identical positive sign (data not shown)]. For the ERS biplot, there is a major separation in two sectors between MN on one side and all other six environments in the other side, with a mean angle vector roughly about a right angle. This biplot view accurately reflects the absence of any significant correlation between MN and the other six locations (Table 8). Regarding the two-by-two relationships among these other six environments as visualized by couples of vectors showing more or less acute angles, they all logically exhibit a significant correlation except three particular couples (LM–SP, GL–SP, and LM–VB). These particular couples of environments had coordinates of opposite sign (negative for SP and VB, positive for GL and LM) on the PC3 axis (data not shown), and the fact that the weight of this axis (14.92%) was marginally less than that of PC2 (15.75%) logically explains the genuine angle values higher than in the two-dimension biplot and therefore loose associations. In the FIB biplot, the percentage of the total GGE variation explained by the first component axis was high (85.52%) as well as the abscissa value of all environment vectors. The combination of these two features suggests coefficients of correlation that are not only significant but also probably high and positive in all cases. This picture was perfectly supported by the correlation values ranging between 0.72 and 0.93 (Table 9). Finally, the EI biplot revealed a large fan-shaped distribution of environment vectors from negative (SB) to positive (LG) ordinates onto a PC2 axis that showed a significant weight (31.23%) hardly inferior to that of PC1 (32.18%). This picture suggests a majority of modest correlation coefficients of either positive or negative sign with only a few, likely positive significant correlations. The values perfectly supported these suggestions (Table 9). Only four positive correlation values were indeed significant (GL–ES, SP–MN, VB–SB, and GL–LM). They logically corresponded for three of them to couples of environments represented by vectors displayed in immediate neighbor positions.

Table 8 Phenotypic correlations of TCH (below diagonal) and ERS (above diagonal) between locations inferred from the balanced dataset used in the GGE study
Table 9 Phenotypic correlations of FIB (below diagonal) and EI (above diagonal) between locations inferred from the balanced dataset used in the GGE study

Mean yield and stability of genotypes

Visualization of both mean performance and stability of genotypes across all environments is always an important issue in cultivar evaluation and to assist in decision-making for selection purposes. This can be done by following the methodology proposed by Yan (2002), which consists in the introduction of an average environment (AE) vector in trait biplots whose coordinates are defined by the average of PC1 and PC2 scores of all environments. This AE vector and its orthogonal vector define a new orthogonal reference frame in which: (i) the highest positive (or negative) abscissa in the new frame should pinpoint the genotypes exhibiting the highest (or lowest) mean performance across all environments, and (ii) the highest ordinates in absolute value (either positive or negative) in this new frame should pinpoint the genotypes exhibiting the highest instability of their performance across environments. These guidelines for analysis in this new reference frame (pictured by dashed arrows in Fig. 2) permits one to identify at a glance: (i) three (G9, G10, G19), one (G2), two (G12, G20), and two (G10, G19) genotypes likely having the highest overall performance means for TCH, ERS, FIB, and EI, respectively, and (ii) one (G3), one (G9), one (G6), and three (G3, G7, G11) genotypes likely having the lowest overall performance means for the same respective traits. All these graphical findings provided by the four two-dimension biplots are perfectly supported by the genotype mean values observed in the original data (Tables 5, 6). They support the reliability of the AE reference frame in displaying a representative ranking of best and poorest mean performances, although biplots did not explain 100% of GGE variation. Moreover, genotypes G11, G7, G7, and G11 appeared to have the highest ordinates in absolute values in this new reference frame for TCH, ERS, FIB, and EI, respectively. These graphical findings pinpoint genotypes that should exhibit the poorest or almost the poorest performance stability across environments. Comparison of these findings with ecovalence values of genotypes (W m ) that measure the contribution of each genotype to the GE interaction (Tables 5, 6) were very congruent. For three of the four traits (TCH, FIB, and EI), the graphical view indeed identified the most unstable genotype, and for the remaining one (ERS) the second most unstable genotype.

Discussion

The present study dissected genotype × environment interactions of traits related to yield components of sugarcane that might exist in the context of sugarcane cultivation on Réunion Island. To this end, we analyzed data from four MET series tested recently in the final stage of eRcane’s selection program based on a network of seven locations representative of the main cultivation zones of sugarcane on Réunion. Genotype-by-environment data of four yield-related traits were investigated to study the main characteristics of genotype response in the multilocation selection scheme. Data were analyzed using mixed linear models to estimate variance components and explored in detail using GGE biplots displaying information on both genotypes and environments.

Importance of the different variance components

Results obtained from our second, “global” model provide estimates of variance components on the basis of the analysis of material belonging to four recent selection series representative of the current breeding program of eRcane. This reference analysis revealed a genetic variance (G) higher than each variance component of the interactions (GL, GC, and GLC) for all traits, except for the key EI trait, for which the GL component was substantially greater than G. Concretely, the GL interaction reflects the mean magnitude of changes in genotype rankings across locations; the GC interaction represents the importance of the fluctuations of performance of genotypes across crop years, which is related to their ratooning ability; the GLC interaction reflects the ratooning ability of genotypes as influenced by location. For all traits, the GL interaction was more important than GC, indicating that testing genotypes across locations is more important than testing for ratooning ability. Besides, TCH and EI produced a GLC highly significantly (P < 0.001) larger than GC. This indicates that the ratooning ability of the genotypes is location specific and illustrates the relative complexity of selecting for these traits, which must be assessed at each location. On the contrary, the null value of the GLC component for FIB indicates that the evolution of FIB across crop years tends to be similar from one location to another.

Simulation of broad-sense heritability (H) with increasing numbers of locations (L) and crop years (C) revealed that: (i) accurate assessment of the mean value of FIB across all environments can be obtained with relatively modest experimental effort (L + C ≤ 3), (ii) the mean performance of genotypes across environments for TCH and ERS required greater experimental effort (L + C > 4) than FIB to reach similar accuracy levels, and (iii) for EI, the lower weight of its genetic variance (G) compared with its GL variance determined H values lower than for the other three traits. For all the traits, the reliability of estimates of genotype means (H values) was influenced more markedly by an increase in the number of locations than by an increase in the number of crop years.

Among the four traits, EI is the most important selection criterion for growers and millers since it represents an economic profitability index. Therefore, to rapidly grasp and compare the potential merit of different candidate genotypes in the diverse cultivation context of Réunion Island, it appears more important to place effort first on the number of trial locations, by identifying additional locations likely to increase the representativeness of the basic location network. In case of budget constraints on resources allocated to field experiments, it would be much more advisable to reduce the number of crop cycles to favor the number of test locations.

Efficiency of GGE plots for visualizing MET data

To evaluate the sugarcane MET network of Réunion Island and obtain essential information for some of the most promising candidate genotypes, we tested the methodology of GGE biplot tools (Yan and Tinker 2006) on a balanced dataset from the two most advanced series (S00 and S03). The two first components of standardized GGE biplots explained a particularly high percentage of initial GGE data for FIB (90%), relatively high for TCH (77%) and EI (72%), and moderately high for EI (63%). Two-dimensional GGE biplots of genotypes and locations sufficiently approximated initial data, providing: (i) a convincing congruency between the comparative arrangement of genotypes on biplot planes of two traits linked by significant positive (TCH and EI) or negative (TCH and ERS) correlations, (ii) easy and reliable visualization of the best performing candidates in each environment, (iii) good immediate visualization of the ranking of genotypes in each environment with an accuracy level directly related to the percentage of initial variation explained for each trait, (iv) a succinct summary of interrelationships among environments, in which acute angles between environment vectors in the same sector depicted at least correlations between environments that are significant when GGE variation is approximately represented by biplots (TCH, ERS, EI), or actually high correlations when variation is represented very efficiently (FIB), and (v) a rapid and broad-brush view of the ranking of mean performance and stability of genotypes across all environments within a new orthogonal reference frame based on a vector of a virtual average environment. All these graphical findings inferred from GGE biplots illustrate the value of these tools to interpret and visualize data results at a glance.

Genotypic response across the multilocation selection network on Réunion Island

Analyses of different genotype-by-location datasets using both mixed models and GGE biplots revealed the general characteristics of genotypic response within the multilocation selection network of Réunion Island and provided valuable information about the selection scheme. Variability between candidate genotypes at this multilocation stage is relatively important for TCH (CV = 20.1%) considering their diverse selection origins from contrasted agroclimatic environments and the original large variation of the trait in earlier selection stages (Dumont pers. commun.). On the contrary, the variability appeared more modest for the ERS (CV = 7.8%) and FIB (CV = 9.8%) quality traits because of the more modest genetic variations than for TCH in the eRcane breeding germplasm and early selection pressure put on ERS since the beginning of selection (Hoarau, pers. commun.). For EI, its variation in this final selection stage appeared relatively high (CV = 20.9%) due to its derived computation from ERS and TCH. Genotype response differed markedly for FIB on one side, for TCH and ERS on a second side, and for EI on a third side. For FIB, the very limited weight of all interaction components (GL, GC, and GLC) compared with the genetic component (G) is reflected by considerably high broad-sense heritability. This implies a high positive correlation between all locations (from 0.72 to 0.92), and many environments are perfectly redundant for genotype ranking for this trait. However, for TCH and ERS, significant coefficients of correlation between environments are much less frequent (10 or 12 occurrences among 21, respectively). Moreover, correlations reached much lower value ranges than for FIB, in alignment with lower heritability levels. Therefore, the chance and strength of redundancy between environments for genotype ranking for each of these two traits are further reduced. Indeed, simultaneous significant correlated couples of environments for both traits (ES/GL, ES/SP, LM/SB, SB/VB) exhibited correlation coefficients in a range not high enough (0.52–0.70) to provide identical or similar genotype ranking among couples of environments for both traits, particularly in the top elites.

Selection strategy for the economic index

EI is the most important criterion for selection for both sugarcane producers and industry, since it represents an economic profitability index. Compared with the other traits, EI is more poorly heritable at the whole MET level (0.70) because of interaction variance components of higher importance. The large part of GGE variation left unexplained by its GGE biplot (more than one-third) combined with the large fan-shaped distribution of environment vectors reflects a nonsignificant or generally loose association between any two environments, as supported by the correlation coefficients (Table 9). This situation indicates an absence of any remarkable similarity of genotype rankings between environments and therefore underlines the importance of the whole selection network. In this context, identification of genotypes adapted to specific environments rather than genotypes with broad adaptation appears to be the most efficient selection strategy to pursue to achieve global genetic progress at the level of the whole sugarcane industry. However, three couples of environments, namely ES/GL, MN/SP, and SB/VB, exhibited associations that were not complete loose, since their correlations for the EI trait ranged from 0.58 to 0.62. Correlations between ES and GL can be explained by the fact that these two sites are separated by a short distance and share the same soil type, differing only in their irrigation regime (Table 1). Correlation between MN and SP can be explained by these sites being relatively unfavorable for cane cultivation with similar chemical fertility of their soil (low pH, low cation exchange capacity) as well as their physical characteristics (very stony), which could imply some similarities in the selection pressure exerted on the tested material. Finally, correlation between SB and VB cannot be obviously explained, since these two environments did not share any similar agroclimatic parameters.

These three particular correlations represent useful information likely to help in rationalization of varieties to be tested in this final MET selection stage. Each year, 5–10 varieties are routinely found to be superior to local controls at each of the seven stations in the previous selection stage. All these elite varieties cannot be tested everywhere due to constraints on land resources and budget. The choice of material entering the MET testing stage always implies difficult decisions between multiple potentially elite genotypes. These modest but interesting correlations between ES and GL, MN and SP, and SB and VB suggest that better genetic gains per unit cost (in the subsequent variety release) could be obtained if greater attention is paid to the choice of a superior proportion of common genotypes between these three couples of environments.

Conclusions

Selection of elite sugarcane genotypes for cultivation on Réunion Island relies on a MET selection program conducted by eRcane research institute. This MET program has been progressively developed over two decades at an increasing number of locations that are representative of the different agroclimatic zones used for sugarcane production. Nowadays, this work routinely includes an unprecedented number of seven experimental stations. Since 2011, each annual new elite series has been systematically tested at all stations of this basic network using a robust standard trial design (in terms of plot size and number of replications) and during an adequate number of crop years to test ratooning ability. Depending on opportunities for collaboration with some producers, some sugarcane series may occasionally be partially tested in a few additional trials at varying locations.

This work is the first report on the MET selection program on Réunion Island in its present-day dimension. We analyzed data so far available from the first four series been tested across the full seven-station network. The choice of the location of these seven stations was guided by the search for the most contrasting areas as possible on the agroclimatic map of the local sugarcane industry, in the hope of efficiently selecting for local adaptation. Data analyses revealed highly significant genotype × location (GL) interaction for all traits of economic interest, as well as for genotype × crop year (GC) interaction, except for FIB (which is determined by a major genetic component). This result confirms a posteriori the interest and relevance of the sites chosen to develop the decentralized selection program likely to exploit or minimize genotype × environment interactions. The optimum number of locations likely to sufficiently assess the mean merit over the whole industry of any elite candidate was determined by plotting the trend of broad-sense heritability (H) for each trait for an increasing number of locations. If sugarcane selection were only based on either tonnes of cane per hectare (TCH) or estimable recoverable sugar (ERS), the marginal increase of H for both of these traits beyond four stations might suggest that not all of them would be necessary. In this hypothetical scenario, resources saved by a reduction in the number of locations could be theoretically guided by considerations relative to the selection of the stations contributing the most to genotype × environment interactions (highest ecovalence index) and the elimination of the very few locations most correlated to one of the former. However, the economic index (EI) is the most important trait for selection (along with resistance to diseases), and its evaluation is based on measurement of both TCH and ERS, two traits that are negatively correlated as observed by other authors (Kang et al. 1983; Milligan et al. 1990, 1996; Baffa et al. 2014). As a result, for this key trait, genotypic response appeared either not correlated between locations or slightly correlated between a few of them. No location appeared redundant relative to the others for this EI criterion, in particular for ranking of top elites. This finding confirms and highlights how selection for local adaptation is an important objective, being more desirable than selection for broad adaptation in the context of sugarcane production on Réunion Island. Moreover, the fact that the simulated broad-sense heritability of EI was not firmly capped by a plateau when reaching seven locations fully justifies the objective of seeking additional opportunities for trials to further test the most promising candidate genotypes (if not their complete series) at supplementary locations.

Our study also illustrated that the first plane of GGE biplots can provide reliable summary representations of responses of candidate genotypes across environments and interrelationships among them. This is directly related to the efficiency of the representation of the variation of the initial trait data. The scope of the lessons that can be drawn from this statistical tool needs to be examined for each trait. In our case study, GGE biplots provide at a glance simple and relatively reliable visualizations of both the performance and stability of genotypes across our MET network when compared with the analytical results. Our case study supports the usefulness of this graphical statistical tool to interpret data results and assist decision-making for selection, further motivating its routine use in sugarcane selection programs.