Adaptability and stability of cotton cultivars (Gossypium hirsutum L. race latifolium H.) using factor analytic model

Abstract

To build up efficient strategies in plant breeding programs, it is requested a certain level of knowledge about the genotype–by-environments interaction (GEI) effects over the crop to be improved. One efficient way to gather this information is using linear mixed models using a parsimonious structure of GEI pattern such as factor analytic (FA) structure. In this work, we applied a multivariate analysis using the FA structure on a dataset composed of 11 cotton genotypes (Gossypium hirsutum) which were evaluated in seven environments in Mozambique, to identify stable and wide/specific adapted genotypes. Using the FA structure, it has been possible to select genotypes with specific and broad adaptability for the environmental network. The results indicated that FK37, Flash, BA525 and BA919 were the most productive genotypes, reaching the highest scores for first Factor. Nevertheless, the FK37 genotype presented Factor 2 close to zero classifying it as the most stable followed by Flash cultivar. The genotypes ISA205, QM301 and Albar SZ9314 were considered the less productive, although QM301 presented good stability. These findings suggest that the FK37 genotype might be recommended for Mozambique, since it had good yield and potential stability when compared to the most cultivated varieties, CA324 and Albar SZ9314.

Introduction

The cotton crop (Gossypium hirsutum L. race latifolium H) is leading in the production of fibres and it is produced commercially worldwide, from temperate to tropical regions (Park et al., 2005; Naveed et al., 2007; Khadi et al., 2010). This crop is fifth in oil production and second in source of protein (Beltrão et al., 2010).

In general, 1 kg of cotton fiber may be obtained starting from 1.65 kg of seed, which contain roughly 21% oil and 23% protein (Wallace et al., 2008; Benzouba et al., 2010). The cotton occupies approximately 34 million hectares distributed across 60 countries around the world, in which Australia (2,000 Kg.ha^-1), Brazil (1,338 Kg.ha^-1), China (1,265 Kg.ha^-1), India (550 Kg.ha^-1), United States of America (985 Kg.ha^-1) and Pakistan (599 Kg.ha^-1) are the most important producers (Fengguo et al., 2007; Beltrão and Azevedo, 2008; Khadi et al., 2010; Acquaah, 2012). The Gossypium genus incorporates 50 different species dispersed worldwide and five of them are tetraploid and belong to the subgenus viz. Karpas.Gossypium hirsutum species is the most cultivated, representing about 96% of cotton world production (Brubaker and Wendel, 1994; Wendel and Cronn, 2003; Fang et al., 2013; Tyagi et al., 2014).

In Mozambique, cotton is one of the most important crops for national economy and contributes about 50 million dollars per year, being classified in the first two most important agricultural exports. Hence, it is a key source of income for more than 300,000 farmers in this country, mainly for small farmers which grades account 90% of the total Mozambique cotton production. In addition, since it is a crop for manufacture, cotton has in its production chain several sectors that employ and/or creating jobs from the farms until the ginning industry (Bias and Donavan, 2003).

One of main challenges of the cotton farmers is the yield variation caused by changes environmental conditions such as related to climate and soils quality that can interfere in crop performance. These factors may cause low heritability or genotypic repeatability which is very common in such quantitative traits. In this scenario, the expected marginal means obtained across several environmental conditions are requested in order to drop out the environmental nuisance factors. Therefore, the evaluation of genotypes in multi-environment trials (METs) network is very important in breeding programs for studying the stability and adaptability of genotypes, as well as its performance prediction in different environments.

In general, GEI may be understood as the phenomenon where the genotypes present different responses across the environments-this scenario may alter the ranking of genotypes in the target environments. Therefore, it is an important issue for breeders and agronomists when the aim is to recommend cultivars for different locations and years or selecting for specific environments or for broad adaptability in mega-environments (Mortazavian et al., 2014). Several methods have been introduced to predict cultivar response in different scenarios (Crossa, 2012; Mortazavian et al., 2014). Among the statistical methods used for MET analysis, we can highlight the linear-bilinear fixed models such us the additive main effects and multiplicative interaction (AMMI) approach and genotype plus GEI (G+GE), the so called GGE biplot (Gauch, 1998; Cornelius and Seyedsadr, 1997; Naveed et al., 2006; Yan et al. 2000). In addition, mixed models approach based on factor analysis or FA structure has been used in MET analysis where the genotypes and GEI are taken as random effects, which genotype and interaction are random effects and environment is fixed (Smith et al., 2015; Stefanova and Buirchell, 2010; Nuvunga et al., 2015). The FA method has been gradually used for analysis in plant breeding (de los Campos and Gianola, 2007; Meyer, 2009; Raman et al., 2011; Zapata-Valenzuela, 2014; Cullis et al., 2014; Smith et al., 2015) and it offers advantages when compared to methods used in the traditional analysis.

Recently, approaches using Bayesian inference for the bilinear models such as AMMI-Bayesian has been applied in MET analysis to overcome some statistical problems related to biplot inferences (Crossa et al., 2011; Josse et al., 2014; Oliveira et al., 2015; Silva et al., 2015). Effective choice of methods for GEI analysis can provide accurate selection of the best varieties for specific or group of environments.

In this study, it was used a factors analysis (FA) to evaluate the response of cotton cultivars (stability and adaptability) in different environments distributed in Mozambique.

Materials and Methods

Experimental details

The experiments were installed in four locations at Namialo (District of Meconta, in Nampula Province), Namara (District of Balama, Cabo Delgado Province), Nhamatanda (District of Nhamatanda, Sofala Province) and Cuamba (District of Cuamba, Niassa Province) during the agricultural seasons 2011/12, 2012/13 and 2013/14, making seven different environments by combining local and years.

In Table 1 relevant characteristics are presented, referring to the different locations that were considered in this analysis.

Table 1. Description agroecologic of locals in evaluation of cotton cultivars productivity in Mozambique.

The relationship between test locations and crop years resulted in seven different environments as shown in Table 2. Thus, our multi-environment dataset is composed of 11 varieties of genotypes (Flash, FK37, BA2018, BA320, BA919, BA525, QM301, ALBAR SZ9314, CA 324 and ISA 205) available in these environments.

Table 2. Environments code combining locals and years.

The experimental design used in each location was a randomized complete block design with four replications. The plots consisted of five lines of 5.0 m long, with the two borders lines with a spacing of 0.70 m between rows and 0.20 cm between plants. The seeds were placed manually, being 4-10seeds per hole with approximately 4 cm deep. The first thinning was done 15 days after emergence leaving two plants per hole and the second thinning was done21 days after emergence leaving a plant per hole. Weeds were controlled manually when needed. The first pest control was done by spraying acetamiprid insecticide (222 g.lt^-1) in a dosage of 50 ml.ha^-1; then was sprayed Lambda-cyhalothrin (60 g.L^-1) fortnightly from the fourth week after emergence, at a dose of 250 ml.ha^-1. Insecticides were applied with micro-ulva (ULV). In the experiment were evaluated the cotton seed yield.

Linear mixed model

A general linear mixed model was used for the first-step analysis in the data from all seven environments. Considering p environments, b blocks and m genotypes, this dataset was subjected to a joint analysis:

y = Xβ +Zu+e (Equation 1)

in which y (n × 1) denotes the combined vector of data across environments, β (pb × 1) are the fixed effects vector (blocks and environments) and u (mp × 1) vector of genotypic effects within locations, X (n × pb) and Z (n × mp) are the associated design matrices, and e (n × 1) denotes the random error vector (m is number of genotypes, b is number of environments and p is number of blocks). It was assumed that e ~ N(0,R) and u ~ N(0,Σ) , in which R is the diagonal residual matrix and Σ is the variance-covariance matrix of genetic effects (Nuvunga et al., 2015).

The factor analytic structure

The genetic variance (Σ) can be modelled through the FA structure, that considers the random effects of m genotypes in p environments as a linear function of latent variables f_k(m × 1) with coefficients λ_k(p × 1) known as loadings, for k=1,..., t < p, and a specific variance vector δ (mp × 1). Then is a matrix (p × t) with the column containing the p environments delivery to the kth latent factor and Ψ is the diagonal matrix of specific variances for the p environments.

The adopted approach (FA-SREG) was described by Meyer et al. (2009) and Nuvunga et al. (2015) based on the FA structure proposed by Smith et al. (2001, 2005), in which the genotypes effects are confounded with GEI, such as the linear-bilinear model of genotypes main effects and interaction GE (Burgueño et al., 2008; Stefanova and Buirchell, 2010). The factors analysis in the restricted maximum likelihood can be easily implemented with the covariance matrix Σ and diagonal R from step one since the unstructured G and diagonal V present very low dimensions. In this situation, the FA can be used as GEI pattern recognition method instead as dimensionality reduction one. It was assumed that G can be represented by an FA structure, such as ΛΛ•+Ψ , and the BLUPs can be represented by common factors as u = Lf+δ, in which . By implementing this transformation in linear mixed model (1) we have that:

y = Xβ +Z[Lf+δ]+e (Equation 2)

in which f (mt × 1) is the factorial scores vector of missing information (BLUP's), δ (mp × 1) is a specific variance vector, L (mp × mt) is the matrix of factorial loadings, and X and Z are the design matrix. Furthermore, it is assumed that f ~ N(0,I) , δ ~ N(0,Ψ) and e ~ N(0,R) .

Therefore, the reparametrized matrix solution to the mixed model equations, in which W = ZL, can be given as follows:

(Equation 3)

Meyer (2009) and Nuvunga et al. (2015) reported the solutions for the fixed and random effects to be similar as follows:

(Equation 4)

(Equation 5)

(Equation 6)

Adaptability and stability of cotton cultivars (Gossypium hirsutum L. race latifolium H.) using factor analytic model

Loadings rotation

When k >1, the loadings matrix ( Λ ) is not unique. This non-uniqueness requires some restriction through the fitting of the FA(k) models because the variance models are termed as non-identifiable (Smith et al., 2002; Cullis et al., 2010).

To ensure the identifiability in Λ, one set of constraints must be imposed on the factor analytic parameters i.e. Λ must be a lower triangular matrix (Smith et al., 2001; Meyer, 2009; Cullis et al., 2010). This is purely for computational simplicity and these restrictions rarely have any basis or biological interpretation. In the context of MET data, a representation of the main components of the loadings (and thus the scores) is usually more significant (Cullis et al., 2010).

After fitting the model, the required rotation is as follows in which:

Λ* = ΛB(Equation 7)

B is the orthogonal matrix, obtained from the singular value decomposition of Λ, whose columns are the eigenvectors of Λ•Λ, More details can be found in Meyer (2009) and Cullis et al. (2010).

Then, means that the first vector (latent variable) accounts for the largest covariance in the FA(k) model, the second vector is orthogonal to the first and it accounts for the next largest amount and so on. In the following, we assume that the loadings correspond to the rotated versions.

Results

The lower triangular unstructured (UN) genotypic (co)variance matrix and diagonal one for residual are shown in Table 3. Based on the Hartley test, there is enough evidence to reject the hypothesis of homogeneity for genetic variance and environmental variance. It suggests the existence of genetic and residual heterogeneity of variances.

Table 3. Genotypic and residual (in red) (co)variance matrix for cotton dataset evaluated in of cotton seed yield at seven different sites.

It is notably also the presence null co-variances among environments (these co-variances represent the genotypic variance plus the variance of interaction between local peers). The pairs of environments with zero covariance can be interpreted as those that present non-correlated genotypic responses contributing largely for GEI. The residual variance was greater for environments 2 and 5, suggesting the influence of those environments in the magnitude of these variances (Table 3).

The FA-SREG2 model explained the 90% of overall genotypic variation. In the Table 4 are the factor loadings rotated by Varimax method. It is observed that the first factor explained almost 80% of the total variation and the second one, 10%. The specific variances were low for all considered environments. The high values of common variance show that the two factors explained a large percentage of each environment variance.

Table 4. Estimated loadings (correlation scale) for the balanced data fit using the FA-SREG2 model.

The Figure 1 shows the factor scores graph of the adjusted FA-SREG2 model. Similarly, to GGEbiplot, the first axis (Factor 1) provides information about the adaptability (the linear correlation of scores related to factor 1 with the marginal BLUPS was 0.90) and the second axis (Factor 2) describes the contribution of the environment or genotype for GEI.

Genotypes whose scores have values greater than zero for Factor 1 are the most productive, while scores near the biplot origin for Factor 2 indicate genotypes that contribute poorly or do not contribute to the GEI; in other words, they are not explained by the environmental loadings and therefore stable. In this scenario, the FK37 genotype was considered the most adapted and stable and it can be widely recommended for all sites in this study. The genotypes that are most distant from the biplot origin for the second score contributed more for the interaction (or specific response) i.e.: Albar SZ9314, BA919, BA525, ISA205.

Figure 2 shows the biplots related to the scores and factorial loadings. In this figure, it is observed four environmental subsets {4, 5}, {1, 3, 6}, {7} and {2}, with respect to the GEI pattern (Figure 2b).

The joint response of genotypes and environments can be examined when in Figure 1. The simultaneous representation of genotypes and environments facilitates interpretation; however, it should be noted that environmental and genotypic scores are not on the same scale. Genotypic and environmental scores located in the same quadrants represent specific positive combinations or adaptability between genotypes and environments.

Discussion

Regarding to the application and interpretation of FA models, there are two issues of primary interest: the estimate of the overall performance of genotypes in all environments and the evaluation of the genotypes stability and adaptability (Stefanova and Buirchell, 2010).

The FA-SREG model can be understood as the linear regression of genotype and GEI in latent environmental covariates (environmental loading up) in which each genotype has a separate slope (genotypic scores, f), but a common intercept (main effects of genotypes are indistinguishable from GEI). The genotypes slopes measure the sensitivity of genotypes to hypothetical environmental factors represented by loading of each environment (Burgueno et al., 2008).

In the FA biplot the first score may represent the general adaptability whether it is correlated with the marginal BLUPs, while the second one may be interpreted as those related to GEI, that is, it evaluates the genotypic stability. Thus, positive values for the first axis correspond to genotypes with high yield and genotypes presenting the second axis values near to biplot origin are those presenting good stability.

Results in Table 3 show the existence of heterogeneity of genotypic and residual variances. A direct implication of this fact would be that the use of analysis of variance (ANOVA)based on least square for GEI and multivariate methods such as AMMI and GGEbiplot (Gauch et al., 2006; Yan et al., 2000; Yan et al., 2007), would not be recommended due the simplification of GEI structure. This claim can also be extended to mixed models or Bayesian methods that simplify the GEI structure using the homogeneity of variance assumption.

Maleia et al. (2010), for instance, evaluated the stability and adaptability of cotton cultivars in Mozambique employing methodologies of Annicchiarico (1992) and of Burrows and Toler (1998). Exploiting the Toler and Burrows (1998) methodology, they identified cotton cultivars presenting specific adaptability for high and low quality environments and others showing wide adaptability for the trail network. In other hand, using the Annicchiarico method, they found cultivars that presented high phenotypic stability.

The genotypes evaluated here were analyzed by Maleia et al. (2017) using fixed AMMI methods found the FK 37, BA 919 e Flash as the most productive genotype and they identified BA2018 e BA320 as more stable genotypes. Although the AMMI method is recognized for allowing the simultaneous study of stability and adaptability, in a single approach, separating pattern of noise, it has the inherent limitations of fixed effects models, such as difficulty in dealing with heteroscedastic data set.

Moiana et al. (2014), on the other hand, using REML/BLUP methodologies (mixed models) to evaluate adaptability and harmonic means to evaluate stability founds genotypes presenting stability and adaptability for Mozambique. Interestingly, the methodology used by these authors is for ranking of the predictions of genotypes and do not evaluate separately the stability and adaptability since the harmonic means of genotypic values is correlated with marginal BLUPs.

However, these methodologies assume the homogeneity of variances between the tested environments which it not requested in mixed models with FA structure. Burgueño et al. (2008) and Cullis et al. (2010) recognized the importance of separating the two GEI causes: the heterogeneity of genetic variances between environments (i.e., the interaction due to scale) and heterogeneity of correlation between environments (Table 3). The last one is often considered more important to MET in plant breeding, since it has impact on the genotypic classification and the selection. The method used in this study allows detailed analysis in this direction, as well as several advantages when compared to traditional methods of analysis as the AMMI fixed effects model and even conventional mixed models such as those used in Maleia et al. (2010), Moiana et al. (2014) and Naveed et al. (2006).

The model used in this study was explained and interpreted the relation of FA-SREG2 and SREG2 models (Burgueno et al. 2008. The biplot interpretation of the of FA models is very similar to the GGEbiplot model presented in Yan et al. (2007), Crossa et al. (2010) and Stefanova and Buirchell (2010). In GGEBiplot, genotypes are evaluated for adaptability from the first principal component scores and in the FA-SREG model this relationship can be obtained through the first factorial score because negative loadings of environments are not observed. In addition, and there was a high correlation between the first factor scores for genotypes and its marginal BLUPs (0.90). The stability of a genotype on GGE analysis, in turn, can be described through the second principal component. A similar interpretation can be given to factor score 2 in the FA-SREG analysis. Smith et al. (2002), in other hand) explained the relation between the FA-AMMI2 and AMMI2 models when the genotype are marginalized from GEI matrix.

In this study, from 11 evaluated genotypes, four are checks- CA324, Albar SZ9314, Chureza and ISA205– being the first two the most grown in Mozambique and the ISA205 is no longer produced. The remaining seven genotypes are genotypes candidates and, therefore, deserve special attention in the trials. The results in Figure 1 of the factorial scores for the FA-SREG2 model shows that the genotypes FK37, Flash, BA525 and BA919 are the most productive reaching the highest scores for Factor 1. Nevertheless, the FK37 genotype with Factor 2 close to zero was the most stable followed by Flash genotype. The genotypes ISA 205, QM301 and Albar SZ9314 are the less productive, although QM301 presents good stability.

As already mentioned, Maleia et al. (2017) analyzed these same genotypes from the conventional AMMI model. In the analysis of adaptability, the results agree with those obtained from the FA analysis. However, with regard to the interaction study, the results diverged between the two approaches. This discrepancy may be partly explained by the fact that we analyzed the MET set in seven environments, while Maleia et al. (2017) evaluated only six of them. In addition, we assume different variances between the test environments, while the AMMI approach has the homogeneity of variances as fundamental assumption.

Conclusion

These results suggest that the FK37 genotype should be recommended for Mozambique, since has a good yield potential and good stability when compared with CA324 and Albar SZ 9314. The FA-SREG model permitted selecting genotypes with specific and broad adaptation in Mozambique.