Integrating Environmental Covariates into Adaptability and Stability Analyses: A Structural Equation Modeling Approach for Cotton Breeding

Abstract

Breeding programs rely on genotype-by-environment interaction (GEI) to recommend cultivars for specific locations. GEI describes how different genotypes perform under varying environmental conditions. Several methods were proposed to assess adaptability and stability across environments. These methods utilize various statistical approaches like parametric and non-parametric regression, multivariate analysis techniques, and even Bayesian frameworks and artificial intelligence. The accessibility of environmental data through platforms like NASA POWER allows breeders to integrate this information into a breeding process. It has been done by using multi-omics integration models that combine data across various biological levels to create accurate predictive models. In the context of phenotypic adaptability and stability analysis, structural equation modeling (SEM) offers an interesting approach to integrating environmental covariates. This work aimed to propose a novel approach that integrates weather information into adaptability and stability analysis, combining SEM with the established Eberhart and Russell model. Additionally, a user-friendly applet, denoted ECERSEM-AdaptStab, was made available to perform the analysis. This approach utilized data from 12 cotton cultivar trials conducted across two growing seasons at 19 sites. This approach successfully integrated environmental covariates into a phenotypic adaptability and stability analysis of cotton cultivars. Specifically, the genotypes TMG 41 WS, IMA CV 690, DP 555 BGRR, BRS 286 and BRS 369 RF were recommended for favorable environments, while the genotypes TMG 43 WS, IMA 5675 B2RF, IMA 08 WS, NUOPAL, DELTA OPAL, BRS 335, and BRS 368 RF are more suitable for unfavorable environments.

1. Introduction

In plant breeding programs, a critical component for recommending cultivars for specific locations is the interaction between the genotypes and the environment (GEI) [1]. GEI refers to the phenomenon where different genotypes perform differently depending on the environmental conditions [2]. Given this factor, particular attention has been given to cotton (Gossypium hirsutum L.), since it is widespread from the southern to the northern hemisphere, from subtropical regions to temperate latitudes [3], with 90% of the total planted area concentrated in the Brazilian Cerrado, a region that is quite large and has 12 Brazilian states [4]. Due to planting in different regions, genotype–environment interactions (GEI) are the main impediment related to the development of cultivars that adapt to environmental heterogeneity and technological specificities [5].

To identify and recommend stable and adapted genotypes across environments, several methodologies to study adaptability and stability have been developed. These methodologies differ in how the adaptability and stability parameters are defined and the methods used for analysis. Some methods rely on parametric and non-parametric regression models, such as those proposed by Finlay and Wilkinson [6] and Eberhart and Russell [7]. Multivariate analysis techniques, including GGE Biplot [8], AMMI [9], and the Extended Centroid Method [10]. Advancements have also seen the Bayesian framework [2,11,12] and artificial intelligence [13] into adaptability and stability analysis.

The advent of platforms like NASA POWER [14] has allowed access to environmental data, making it easier to collect this information and use it than in breeding programs. Specifically, it is allowing breeders to integrate environmental covariates and genomic information through a multi-omics approach [15,16].

A natural way to integrate environmental covariates into adaptability and stability analyses is through structural equation modeling (SEM). SEM allows to estimate both direct and indirect effects within a system of interrelated variables, assuming a linear model. This powerful capability paves the way for integrating weather data into adaptability and stability analyses. By combining SEM with the established Eberhart and Russell model, we can incorporate weather covariates. Specifically, this approach can allow to estimate the direct impact of weather (climate covariates) on the environmental index used in the Eberhart and Russell model, as well as its direct and indirect influence on the yield of the evaluated genotypes.

SEM has been used successfully in breeding to infer causal phenotype networks using structural equation models [17], investigate multi-trait genetic architecture of udder health in dairy cattle [18], and in GWAS studies to explore interrelated dependencies between phenotypes related to morphology in coffee [19]. However, to the best of our knowledge, this is the first attempt to leverage SEM within the adaptability and stability framework for crop breeding.

In light of this, the purposes of this study were to: (i) propose a novel approach that integrates weather information into combing structural equation modeling (SEM) with the established Eberhart and Russell model; (ii) make available an applet (Adaptability and Stability Analysis incorporating Environmental Covariates into Eberhart and Russell based on Structural Equation Models—ECERSEM-AdaptStab) that allows us to collect the weather information and provide for the user both direct and indirect effects from the climate covariates in the environmental index and on the yield of the evaluated genotypes.

2. Materials and Methods

2.1. Experimental Data Set

The data set was obtained from a cotton cultivar competition trial with 12 treatments planted in 19 environments following a randomized block design with four replications in the years 2013/2014 and 2014/2015 (Figure 1). The 12 treatments were composed of the following mid-cycle cultivars: TMG 41 WS (genotype 1), TMG 43 WS (genotype 2), IMA CV 690 (genotype 3), IMA 5675 B2RF (genotype 4), IMA 08 WS (genotype 5), NUOPAL (genotype 6), DP 555 BGRR (genotype 7), DELTA OPAL (genotype 8), BRS 286 (genotype 9), BRS 335 (genotype 10), BRS 368 RF (genotype 11), and BRS 369 RF (genotype 12). The environments were produced from the combination of 2 years and different production locations in the Brazilian Cerrado, more specifically the regions of Minas Gerais (MG), Goiás (GO), Mato Grosso (MT), Bahia (BA), Maranhão (MA), and Piauí (PI) (Figure 1). Due to the vast size of this territory, it presents different edaphoclimatic characteristics, which can be justified by the latitudinal variation, topographic amplitude, inland position, and the action of air masses [20,21,22]. All experiments had an experimental unit composed of four 5 m lines with spacings of 0.90 m between the lines and a density of 100,000 $\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{t}\mathrm{s}.{\mathrm{h}\mathrm{a}}^{-1}$. The trait evaluated was fiber yield. The collection was carried out from samples from the two central rows of each experimental unit, and finally the humidity was corrected to 13%. More details about the experiment can be found in Teodoro et al. [23].

In order to understand the effect of environmental covariates on the performance of each genotype in each environment, the covariates related to Earth’s skin temperature (TS), temperature at 2 meters (T2M), relative humidity at 2 meters (RH2M), wind speed at 2 meters maximum (WS2M_MAX), corrected precipitation (PRECTOTCORR), and total clear sky surface photosynthetically active radiation (CLRSKY_SFC_PAR_TOT) were obtained from the NASA POWER platform (https://power.larc.nasa.gov/data-access-viewer/, accessed on 15 September 2024) [14].

2.2. Statistical Analysis

Initially, the joint analysis variance (ANOVA) was performed to verify if the effect of GEI is significant according to the model:

$${\mathrm{Y}}_{\mathrm{i}\mathrm{j}\mathrm{k}}=\mathsf{\mu }+{\mathrm{G}}_{\mathrm{i}}+{\mathrm{E}}_{\mathrm{j}}+{\mathrm{B}/\mathrm{E}}_{\mathrm{j}\mathrm{k}}+{\mathrm{G}\mathrm{E}\mathrm{I}}_{\mathrm{i}\mathrm{j}}+{\mathrm{e}}_{\mathrm{i}\mathrm{j}\mathrm{k}},$$

(1)

where ${\mathrm{Y}}_{\mathrm{i}\mathrm{j}\mathrm{k}}$ is the observation in the $\mathrm{k}$th block, evaluated in the $\mathrm{j}$th environment and $\mathrm{i}$th genotype; $\mathsf{\mu }$ is the overall mean; ${\mathrm{G}}_{\mathrm{i}}$ is the effect of the $\mathrm{i}$th genotype (fixed effect); ${\mathrm{E}}_{\mathrm{j}}$ is the effect of the $\mathrm{j}$th environment (fixed effect); ${\mathrm{B}/\mathrm{E}}_{\mathrm{j}\mathrm{k}}$ is the effect of block $\mathrm{k}$ within environment $\mathrm{j}$ (fixed effect); ${\mathrm{G}\mathrm{E}\mathrm{I}}_{\mathrm{i}\mathrm{j}}$ is the effect of the genotype $\mathrm{i}$ x environment $\mathrm{j}$ interaction (fixed effect); and ${\mathrm{e}}_{\mathrm{i}\mathrm{j}\mathrm{k}}$ is the error associated with the ${\mathrm{Y}}_{\mathrm{i}\mathrm{j}\mathrm{k}}$ observation (random effect). All this analysis was carried out using the R software version 4.3.1 [24].

The traditional Eberhart and Russell methodology [7] is based on simple linear regression:

$${\mathrm{Y}}_{\mathrm{i}\mathrm{j}}={\mathsf{\beta }}_{0\mathrm{i}}^{\mathrm{E}\mathrm{R}}+{\mathsf{\beta }}_{1\mathrm{i}}^{\mathrm{E}\mathrm{R}}{\mathrm{I}}_{\mathrm{j}}+{\mathsf{\delta }}_{\mathrm{i}\mathrm{j}}+{\mathrm{e}}_{\mathrm{i}\mathrm{j}},$$

(2)

where $\mathrm{Y}$_ij is the observation of the $\mathrm{i}$th (i = 1, 2, …, g) genotypes in the $\mathrm{j}$th (j = 1, 2, …n) environment; ${\mathsf{\beta }}_{0\mathrm{i}}^{\mathrm{E}\mathrm{R}}$ is the general mean for the $\mathrm{i}$th genotype; ${\mathsf{\beta }}_{1\mathrm{i}}^{\mathrm{E}\mathrm{R}}$ is the regression coefficient for the $\mathrm{i}$th genotype; $\mathrm{I}$_j is the environmental index $\left({\mathrm{I}}_{\mathrm{j}}={\displaystyle \frac{\sum _{\mathrm{i}}{\mathrm{Y}}_{\mathrm{i}\mathrm{j}}}{\mathrm{g}}}-{\displaystyle \frac{\sum _{\mathrm{i}}\sum _{\mathrm{j}}{\mathrm{Y}}_{\mathrm{i}\mathrm{j}}}{\mathrm{g}\mathrm{a}}}\right)$; ${\mathsf{\delta }}_{\mathrm{i}\mathrm{j}}$ is equal to the regression deviation of the $\mathrm{i}$th genotype in the $\mathrm{j}$th environment, and $\mathrm{e}$_ij is the effect of the mean experimental error. In this approach, genotype performance refers to how each genotype performs across different environmental conditions, and it is measure by the regression parameter, ${\mathsf{\beta }}_{1\mathrm{i}}^{\mathrm{E}\mathrm{R}}$. Stability parameter, ${\mathsf{\sigma }}_{\mathrm{d}\mathrm{S}\mathrm{E}\mathrm{M}}^2$, on the other hand, informs us about the predictability of a genotype’s behavior in response to environmental stimuli. The stability can be measure also considering the coefficient of determination [25].

2.3. Structural Equation Models

Aiming to integrate weather data into adaptability and stability analyses, a structural equation modeling framework was proposed. The causal relationship between the endogenous variables (average yield of each genotype across the environments and the environmental index—$\mathbf{I}$), as well as the exogenous variables (environmental covariables), was modeled as follows [26].

$${\mathbf{y}}^{*}=\mathsf{\alpha }+\mathbf{B}{\mathbf{y}}^{*}+\mathsf{\Gamma }\mathbf{x}+\mathsf{\zeta },$$

(3)

where ${\mathbf{y}}^{*}$ (average fiber yield of each genotype across the environments and the environmental index) and $\mathbf{x}$ (TS, T2M, RH2M, WS2M_MAX, PRECTOTCORR and CLRSKY_SFC_PAR_TOT) are (2 × 1) and (2 × 6) vectors that represent the endogenous and exogenous variables, respectively. $\mathsf{\alpha }$ is a (2 × 1) vector that represents the general mean of each genotype (general mean, as Eberhart and Russell). $\mathbf{B}$ and $\mathsf{\Gamma }$ are (2×2) and (6×1) matrices that represent the causal coefficients representing direct effects of $\mathbf{y}$-variables on ${\mathbf{y}}^{*}$-variables (adaptability parameter, as Eberhart and Russell) and direct effects of $\mathbf{x}$-variables on ${\mathbf{y}}^{*}$-variables, respectively, and $\mathsf{\zeta }$ is a (2 × 1) vector that represents the random residuals associated with endogenous variables. Using the SEM theory, it was possible to verify indirect effects of EC via the environmental index in each of the genotypes. These values are defined by the product between the environmental covariable estimates and the adaptability parameter (${\widehat{\mathsf{\gamma }}}_{\mathrm{j}}\times {\widehat{\mathsf{\beta }}}_{1\mathrm{i}}^{\mathrm{S}\mathrm{E}\mathrm{M}}$). The estimation for SEM model was based on maximum likelihood estimation (ML). According to Boolen [27], the coefficient of determination is given by ${\mathrm{R}}_{\mathrm{S}\mathrm{E}\mathrm{M}}^2=1-{\displaystyle \frac{\left|\widehat{\mathsf{\Psi }}\right|}{\left|{\widehat{\mathsf{\Sigma }}}_{\mathbf{y}\mathbf{y}}\right|}}$, where $\left|\widehat{\mathsf{\Psi }}\right|$ is the determinant of the estimate of the covariance matrix for the equation errors and $\left|{\widehat{\mathsf{\Sigma }}}_{\mathbf{y}\mathbf{y}}\right|$ is the determinant of the covariance matrix of ${\mathbf{y}}^{*}$ taken from $\widehat{\mathsf{\Sigma }}$ (estimated covariance matrix). The stability parameter is estimated by ${\widehat{\mathsf{\sigma }}}_{\mathrm{d}\mathrm{S}\mathrm{E}\mathrm{M}}^2={\mathsf{\psi }}_{11}=\mathrm{v}\mathrm{a}\mathrm{r}\left({\mathrm{y}}_{\mathrm{i}}\right)-{\mathsf{\sigma }}_{{\mathrm{y}}_{\mathrm{i}}{\mathrm{z}}_{\mathrm{i}}}^{\prime }{\mathsf{\Sigma }}_{{\mathrm{z}}_{\mathrm{i}}{\mathrm{z}}_{\mathrm{i}}}^{-1}{\mathsf{\sigma }}_{{\mathrm{z}}_{\mathrm{i}}{\mathrm{y}}_{\mathrm{i}}}$, where $\mathrm{v}\mathrm{a}\mathrm{r}\left({\mathrm{y}}_{\mathrm{i}}\right)$ is the variance of ${\mathrm{y}}_{\mathrm{i}}$, ${\mathsf{\Sigma }}_{{\mathrm{z}}_{\mathrm{i}}{\mathrm{z}}_{\mathrm{i}}}^{-1}$ is the nonsingular covariance matrix of the column vector ${\mathrm{z}}_{\mathrm{i}}$ contains only those variables that belong in the ${\mathrm{y}}_{\mathrm{i}}$ equation, and ${\mathsf{\sigma }}_{{\mathrm{z}}_{\mathrm{i}}{\mathrm{y}}_{\mathrm{i}}}$ is the covariance vector of ${\mathsf{\zeta }}_1$ (random residuals associated with ${\mathrm{y}}_{\mathrm{i}}$). These analyses were performed using the lavaan package [28] in R software [24].

The causal relationship considering the proposed model can be written individually as:

$$\begin{gathered} \mathbf{I}={\mathsf{\gamma }}_{\mathrm{T}\mathrm{S}\to \mathrm{I}}\mathrm{T}\mathrm{S}+{\mathsf{\gamma }}_{\mathrm{T}2\mathrm{M}\to \mathrm{I}}\mathrm{T}2\mathrm{M}+{\mathsf{\gamma }}_{\mathrm{R}\mathrm{H}2\mathrm{M}\to \mathrm{I}}\mathrm{R}\mathrm{H}2\mathrm{M}+{\mathsf{\gamma }}_{\mathrm{W}\mathrm{S}2\mathrm{M}\_\mathrm{M}\mathrm{A}\mathrm{X}\to \mathrm{I}}\mathrm{W}\mathrm{S}2\mathrm{M}\_\mathrm{M}\mathrm{A}\mathrm{X}+{\mathsf{\gamma }}_{\mathrm{P}\mathrm{R}\mathrm{E}\mathrm{C}\mathrm{T}\mathrm{O}\mathrm{T}\mathrm{C}\mathrm{O}\mathrm{R}\mathrm{R}\to \mathrm{I}}\mathrm{P}\mathrm{R}\mathrm{E}\mathrm{C}\mathrm{T}\mathrm{O}\mathrm{T}\mathrm{C}\mathrm{O}\mathrm{R}\mathrm{R}+ \\ {\mathsf{\gamma }}_{\mathrm{C}\mathrm{L}\mathrm{R}\mathrm{S}\mathrm{K}\mathrm{Y}\_\mathrm{S}\mathrm{F}\mathrm{C}\_\mathrm{P}\mathrm{A}\mathrm{R}\_\mathrm{T}\mathrm{O}\mathrm{T}\to \mathrm{I}}\mathrm{C}\mathrm{L}\mathrm{R}\mathrm{S}\mathrm{K}\mathrm{Y}\_\mathrm{S}\mathrm{F}\mathrm{C}\_\mathrm{P}\mathrm{A}\mathrm{R}\_\mathrm{T}\mathrm{O}\mathrm{T}+{\mathsf{\zeta }}_2, \end{gathered}$$

(4)

$$\begin{gathered} {\mathrm{y}}_{\mathrm{i}}={\mathsf{\beta }}_{0\mathrm{i}}^{\mathrm{S}\mathrm{E}\mathrm{M}}+{\mathsf{\gamma }}_{\mathrm{T}\mathrm{S}\to {\mathrm{y}}_{\mathrm{i}}}\mathrm{S}\mathrm{T}+{\mathsf{\gamma }}_{\mathrm{T}2\mathrm{M}\to {\mathrm{y}}_{\mathrm{i}}}\mathrm{T}2\mathrm{M}+{\mathsf{\gamma }}_{\mathrm{R}\mathrm{H}2\mathrm{M}\to {\mathrm{y}}_{\mathrm{i}}}\mathrm{R}\mathrm{H}2\mathrm{M}+{\mathsf{\gamma }}_{\mathrm{W}\mathrm{S}2\mathrm{M}\_\mathrm{M}\mathrm{A}\mathrm{X}\to {\mathrm{y}}_{\mathrm{i}}}\mathrm{W}\mathrm{S}2\mathrm{M}\_\mathrm{M}\mathrm{A}\mathrm{X}+ \\ {\mathsf{\gamma }}_{\mathrm{P}\mathrm{R}\mathrm{E}\mathrm{C}\mathrm{T}\mathrm{O}\mathrm{T}\mathrm{C}\mathrm{O}\mathrm{R}\mathrm{R}\to {\mathrm{y}}_{\mathrm{i}}}\mathrm{P}\mathrm{R}\mathrm{E}\mathrm{C}\mathrm{T}\mathrm{O}\mathrm{T}\mathrm{C}\mathrm{O}\mathrm{R}\mathrm{R}+{\mathsf{\gamma }}_{\mathrm{C}\mathrm{L}\mathrm{R}\mathrm{S}\mathrm{K}\mathrm{Y}\_\mathrm{S}\mathrm{F}\mathrm{C}\_\mathrm{P}\mathrm{A}\mathrm{R}\_\mathrm{T}\mathrm{O}\mathrm{T}\to {\mathrm{y}}_{\mathrm{i}}}\mathrm{C}\mathrm{L}\mathrm{R}\mathrm{S}\mathrm{K}\mathrm{Y}\_\mathrm{S}\mathrm{F}\mathrm{C}\_\mathrm{P}\mathrm{A}\mathrm{R}\_\mathrm{T}\mathrm{O}\mathrm{T}+{\mathsf{\beta }}_{1\mathrm{i}}^{\mathrm{S}\mathrm{E}\mathrm{M}}\mathbf{I}+{\mathsf{\zeta }}_1, \end{gathered}$$

(5)

where $\mathbf{I}$ is the vector of the environmental index, ${\mathrm{y}}_{\mathrm{i}}$ (average fiber yield of each genotype across the environments) is the vector of the $\mathrm{i}$th genotype ($\mathrm{i}=1,\dots ,12)$, ${\mathsf{\beta }}_{0\mathrm{i}}^{\mathrm{S}\mathrm{E}\mathrm{M}}$ is the vector that represent the general mean of each genotype (general mean, as Eberhart and Russell). $\mathsf{\gamma }$ values represent the direct effects among environmental covariates (TS, T2M, RH2M, WS2M_MAX, PRECTOTCORR, and CLRSKY_SFC_PAR_TOT) on the index and each genotype. The ${\mathsf{\beta }}_{1\mathrm{i}}^{\mathrm{S}\mathrm{E}\mathrm{M}}$ values represent the direct effect among environmental indexes on each genotype (adaptability parameter, as Eberhart and Russell). ${\mathsf{\zeta }}_1$ and ${\mathsf{\zeta }}_2$ represents the random residuals associated with endogenous variables. Also, the stability can be measured considering the coefficient of determination [25]. The causal relationships represented by Equations (4) and (5) of the structural equation model are visualized in the causal diagram of Figure 2.

2.4. Integrating Weather Data into Adaptability and Stability Analyses into a Shiny Applet

All the aforementioned methodologies were integrated into the Adaptability and Stability Analysis incorporating Environmental Covariates into Eberhart and Russell based on Structural Equation Models (ECERSEM-AdaptStab). The ECERSEM-AdaptStab is a freely available tool at https://ecersem-adaptstab.shinyapps.io/ECERSEM-AdaptStab/ (accessed on 21 October 2024) and offers five functionalities: (1) data upload and covariate selection; (2) ANOVA for GEI significance; (3) data summary; (4) heatmaps for adaptability and stability parameters (${\widehat{\mathsf{\beta }}}_{0\mathrm{i}}^{\mathrm{S}\mathrm{E}\mathrm{M}}$, ${\widehat{\mathsf{\beta }}}_{1\mathrm{i}}^{\mathrm{S}\mathrm{E}\mathrm{M}}$, ${\widehat{\mathsf{\sigma }}}_{\mathrm{d}\mathrm{S}\mathrm{E}\mathrm{M}}^2$ and ${\mathrm{R}}_{\mathrm{S}\mathrm{E}\mathrm{M}}^2$) and; (5) heatmaps for environmental effects on genotypes and index. The workflow for using the ECERSEM-AdaptStab app is presented in Figure 3.

For part 1: (i) to ensure compatibility with the ECERSEM-AdaptStab app, the phenotypic information file must have the following format: columns named “LOCAL”, “TRAT”, “REP”, and “VARIABLE” in that specific order. It is important that the information in the first three columns is numeric and starts from “1” (for example, if the phenotypic database has 3 sites, 5 treatments, and 4 replications, we would have a table where the first site, first treatment, and first replication must be entered as “1”, going up to the values “3”, “5”, and “4”, respectively, for this hypothetical situation). (ii) The local information file should also follow this specific format: five columns named “LOCAL”, “LATITUDE”, “LONGITUDE”, “DATE_INITIAL”, and “DATE_FINAL” in that exact order. Additionally, the “LOCAL” column should contain unique identifiers (usually text labels). The user must use the decimal degrees format for latitude and longitude and “YEAR-MM-DAY” format for dates (Figure 4).

Parts 4 and 5 of the tools present heatmaps summarizing the results of the SEM analysis. Importantly, these analyses are based on a standard SEM structure, which is visualized in Figure 3. It is observed that the model contemplates the influence of environmental covariables (TS, T2M, RH2M, WS2M_MAX, PRECTOTCORR, and CLRSKY_SFC_PAR_TOT) in average fiber yield of each genotype (${\mathrm{y}}_{\mathrm{i}}$) and the environmental index ($\mathbf{I}$) based on the $\widehat{\mathsf{\gamma }}$ values, and the influence of the environmental index on each genotype (adaptability parameter) based on the $\widehat{\mathsf{\beta }}$ values.

3. Results

A summary of the descriptive statistics, including the means, medians, standard deviations (SD), and ranges (minimum and maximum) for the fiber yield of each genotype, the environmental index, and each environmental covariate (EC), is presented in

Table 1. Means, medians, standard deviations (SD), and ranges (minimum and maximum) for the fiber yield of each genotype, the environmental index, and each environmental covariate.

Variables	Min	Max	Mean	Median	SD
FYG1	499.86	2416.36	1697.10	1759.41	508.19
FYG2	801.66	2297.76	1683.81	1698.99	438.19
FYG3	587.90	3093.16	1951.87	1873.70	591.22
FYG4	857.91	2220.70	1690.92	1777.35	355.78
FYG5	791.17	2473.47	1764.82	1830.74	410.74
FYG6	697.39	2585.81	1645.67	1670.49	471.06
FYG7	483.95	2786.23	1899.76	1876.04	599.87
FYG8	572.91	2592.31	1653.66	1691.46	564.53
FYG9	500.61	2604.90	1732.90	1701.33	516.39
FYG10	544.79	2614.17	1688.66	1813.26	562.52
FYG11	808.54	2889.71	1773.73	1805.11	540.02
FYG12	631.40	2712.45	1818.07	1914.97	571.76
Environmental index	−1083.63	735.94	0.00	87.20	469.92
TS	23.15	26.13	24.60	24.29	1.02
T2M	23.06	26.00	24.32	24.06	0.98
RH2M	65.77	90.38	78.81	80.29	6.74
WS2M_MAX	0.16	3.73	2.66	3.18	1.17
PRECTOTCORR	1.28	9.76	5.16	5.42	3.19
CLRSKY_SFC_PAR_TOT	141.16	156.78	151.32	152.30	3.70

FYG: average fiber yield of each genotype (1, 2, …, 12); TS: Earth skin temperature; T2M: temperature at 2 meters; RH2M: relative humidity at 2 meters; WS2M_MAX: wind speed at 2 meters maximum; PRECTOTCORR: precipitation corrected; CLRSKY_SFC_PAR_TOT: clear sky surface photosynthetically active radiation total. Min: minimum; Max: maximum; SD: standard deviation.

. The average of fiber yield values varied from 1645.67 (±471.06) kg$.{\mathrm{h}\mathrm{a}}^{-1}$(NUOPAL) to 1899.76 (±599.87) kg$.{\mathrm{h}\mathrm{a}}^{-1}$ (genotype 7). The minimum and maximum values ranged from 483.95 kg$.{\mathrm{h}\mathrm{a}}^{-1}$ (DP 555 BGRR) to 3093.16 kg$.{\mathrm{h}\mathrm{a}}^{-1}$ (IMA CV 690). Eleven of the nineteen environments (1, 2, 3, 6, 7, 9, 12, 13, 16, 17, and 19) presented positive values and eight (4, 5, 8, 10, 11, 14, 15, and 18) negative values for the environmental index, with values ranging from 735.94 to −1083.63. The mean values of the EC, TS, T2M, RH2M, WS2M_MAX, PRECTOTCORR, and CLRSKY_SFC_PAR_TOT presented values ranging from 24.60 (±1.02) $\mathrm{℃}$, 24.32 (±0.98) $\mathrm{℃}$, 78.81 (±6.74) %, 2.66 (±1.17) ${\mathrm{m}.\mathrm{s}}^{-1}$, 5.16 (±3.19) ${\mathrm{m}\mathrm{m}.\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}$ and 151.32 (±3.70) ${\mathrm{W}.\mathrm{m}}^{-2}$, respectively (Table 1).

3.1. Analysis of Variance for Cotton Yields in Different Environments

The joint analysis of variance for the data of yield of 12 cotton genotypes showed the significance (p-value $<$ 0.01) of the GEI (

Table 2. Joint analysis for yield of 12 early cotton genotypes evaluated in 19 Brazilian environments in the 2013/2014 and 2014/2015 crop seasons.

Source of Variation	Degrees of Freedom	Mean Square
Genotype	11	714,991.03 *
Environments	18	10,599,770.42 *
Block/environments	57	42,223.61 *
GEI	198	198,986.84 *
Residual	627	23,035.85
Mean	-	1750.08
Coefficient of variation (%)	-	8.67

* Represents significance at 1% probability by the F test.

). This result indicates the differential performance of genotypes in different environments, justifying the use of adaptability and stability analysis. Additionally, the average fiber yield was 1750.08 kg$.{\mathrm{h}\mathrm{a}}^{-1}$ and the coefficient of variation obtained was 8.67%.

3.2. Yield Adaptability and Stability from Experimental Data

Observing the ECERSEM-AdaptStab and ER results, we can see that the genotypes IMA CV 690, IMA 08 WS, DP 555 BGRR, BRS 368 RF, and BRS 369 RF presented ${\widehat{\mathsf{\beta }}}_0$ higher than the general mean (1750.08) (

Table 3. Adaptability and stability parameters obtained by the ECERSEM-AdaptStab app for 12 cotton genotypes in 19 environments in the Brazilian Cerrado.

	ECERSEM-AdaptStab				ER
Genotypes	${\widehat{\mathsf{\beta }}}_0^{\mathbf{S}\mathbf{E}\mathbf{M}}$	${\widehat{\mathsf{\beta }}}_1^{\mathbf{S}\mathbf{E}\mathbf{M}}$	${\widehat{\mathsf{\sigma }}}_{\mathbf{d}\mathbf{S}\mathbf{E}\mathbf{M}}^2$	${\mathbf{R}}_{\mathbf{S}\mathbf{E}\mathbf{M}}^2$	${\widehat{\mathsf{\beta }}}_0^{\mathbf{E}\mathbf{R}}$	${\widehat{\mathsf{\beta }}}_1^{\mathbf{E}\mathbf{R}}$	${\widehat{\mathsf{\sigma }}}_{\mathbf{d}\mathbf{E}\mathbf{R}}^2$	${\mathbf{R}}_{\mathbf{E}\mathbf{R}}^2$
TMG 41 WS	1697.10	1.09 *	29,027.54 *	0.88	1697.10	0.98	42,531.17 *	0.82
TMG 43 WS	1683.81	0.91 *	40,204.28 *	0.77	1683.81	0.80 *	47,251.75 *	0.74
IMA CV 690	1951.87	1.52 *	14,313.06 *	0.95	1951.87	1.19 *	30,820.69 *	0.90
IMA 5675 B2RF	1690.92	0.92 *	19,863.71 *	0.83	1690.92	0.67 *	24,638.68 *	0.77
IMA 08 WS	1764.82	0.70 *	39,891.97 *	0.75	1764.82	0.70 *	58,017.47 *	0.64
NUOPAL	1645.67	0.43 *	9471.79 *	0.95	1645.67	0.94	22,563.10 *	0.88
DP 555 BGRR	1899.76	1.55 *	5824.48 *	0.98	1899.76	1.24 *	12,839.75 *	0.95
DELTA OPAL	1653.66	0.92 *	21,167.75 *	0.92	1653.66	1.13 *	35,320.60 *	0.88
BRS 286	1732.90	1.25 *	22,869.41 *	0.91	1732.90	1.03	30,571.81 *	0.87
BRS 335	1688.66	0.52 *	35,714.77 *	0.88	1688.66	1.08 *	58,215.95 *	0.81
BRS 368 RF	1773.73	0.78 *	30,656.38 *	0.88	1773.73	1.06	38,783.87 *	0.86
BRS 369 RF	1818.07	1.36 *	12,435.91 *	0.95	1818.07	1.18 *	15,010.33 *	0.94
General mean	1750.08	-	-	-	1750.08	-	-	-

ECERSEM-AdaptStab: adaptability and stability analysis incorporating environmental covariates into Eberhart and Russell based on structural equation models; ER: Eberhart and Russell; * represents significance at 5% probability by the F test.

). According to ECERSEM-AdaptStab, the genotypes TMG 41 WS, IMA CV 690, DP 555 BGRR, BRS 286, and BRS 369 RF presented ${\mathsf{\beta }}_1>1$ and are recommended for favorable environments. On the other hand, the genotypes TMG 43 WS, IMA 5675 B2RF, IMA 08 WS, NUOPAL, DELTA OPAL, BRS 335, and BRS 368 RF were classified as genotypes with specific adaptability to unfavorable environments (${\mathsf{\beta }}_1<1$). None of the evaluated genotypes were classified as genotypes of general adaptability ${(\mathsf{\beta }}_1=1)$. Among the evaluated genotypes, observing ECERSEM-AdaptStab only TMG 43 WS and IMA 08 WS, and in ER TMG 43 WS, IMA 5675 B2RF, and IMA 08 WS presented low stability, ${\mathrm{R}}^2<0.80$ (Table 3). The results using ER show that IMA CV 690, DP 555 BGRR, DELTA OPAL, BRS 335, and BRS 369 RF presented ${\mathsf{\beta }}_1>1$, TMG 43 WS, IMA 5675 B2RF, and IMA 08 WS presented ${\mathsf{\beta }}_1<1$, and TMG 41 WS, NUOPAL, BRS 286, and BRS 368 RF presented ${\mathsf{\beta }}_1=1$.

The direct effects from EC on genotypes are presented in Figure 5 and in table format in Table S1. Among the evaluated genotypes, TMG 41 WS, TMG 43 WS, IMA 5675 B2RF, IMA 08 WS, and BRS 368 RF did not exhibit significant interrelationships (Figure 5). The environmental covariates TS and T2M showed larger magnitude differences between the negative and positive effects. The magnitude of the influence is indicated by the values in parentheses. TS presented a negative impact on IMA CV 690 (−636.36), DP 555 BGRR (−354.86), and BRS 286 (−550.70), while positively affecting genotypes NUOPAL (541.66), DELTA OPAL (666.4), and BRS 335 (646.6). Similar to TS, T2M negatively influenced genotypes NUOPAL (−455.39) and DELTA OPAL (−622.49) and positively influenced IMA CV 690 (490.07), DP 555 BGRR (389.41), and BRS 286 (538.3). Only three genotypes were significantly influenced by RH2M. IMA CV 690 (−30.01) was negatively affected and DELTA OPAL (42.61) and BRS 369 RF (20.84) responded positively. Wind speed (WS2M_MAX) negatively affected DP 555 BGRR (−118.73) but positively influenced NUOPAL (135.71) and BRS 335 (212.70). PRECTOTCORR presented a positive effect on IMA CV 690 (14.8). Finally, CLRSKY_SFC_PAR_TOT showed a negative influence only on NUOPAL (−25.02) and a positive influence on IMA CV 690 (30.17) and DP 555 BGRR (46.06).

The indirect effects of EC on the genotype via environmental index are presented in Figure 6. TS and WS2M_MAX positively influenced the environmental index with values of 795.90 and 321.30, respectively. T2M and CLRSKY_SFC_PAR_TOT negatively influenced the environmental index with values of −621.90 and −55.30, respectively (Figure 5, Table S1).

According to Figure 7, analyzing TS and WS2M_MAX, all genotypes showed positive total effects, with values ranging from 450.05 to 1403.84 and 225.82 to 489.14, respectively. Evaluating T2M and CLRSKY_SFC_PAR_TOT, it was observed that all genotypes showed negative total effects, with values ranging from −1198.65 to −243.58 and −75.24 to −28.81, respectively. For RH2M and PRECTOTCORR, not all genotypes showed significant effects, as the index did not significantly affect these same ECs (Figure 5). For these ECs, the total effect on genotypes was exactly the same as the direct effect (Figure 5).

4. Discussion

In this study, we proposed a structural equation modeling (SEM) approach to integrate weather data into analyses of genotype adaptability and stability. Beyond estimating adaptability and stability parameters like the Eberhart and Russell approach, SEM allows us to assess both the direct and indirect effects of environmental covariates on genotypes through the environmental index. To present the proposal, we utilized data from a cotton cultivar competition trial with 12 treatments planted across 19 environments. Additionally, a user-friendly ECERSEM-AdaptStab app for performing this analysis was developed and made available.

Similar to several studies performed to investigate cotton crop adaptability and stability in the Brazilian Cerrado, the GEI was significant for the same trait in regions similar to this study [23,25,29,30,31]. The observed coefficient of variation (CV%) for cotton fiber yield (8.67%) is consistent with values reported in the literature [12,32] and is indicating high experimental precision [33].

Initially, the results presented from ECERSEM-AdaptStab were compared with those obtained using the traditional ER methodology (Table 3). This comparison revealed that all genotypes exhibited convergent results for both the overall mean (general average) and the stability parameter based on component of variance (${\widehat{\mathsf{\sigma }}}_{\mathrm{d}}^2$). Regarding the stability parameter measured by the coefficient of determination (${\mathrm{R}}^2$), only the IMA 5675 B2RF genotype showed a discrepancy. ECERSEM-AdaptStab estimated its stability to be above 80%, whereas the traditional ER method yielded a value below 80%. The genotypes TMG 43 WS, IMA CV 690, IMA 5675 B2RF, IMA 08 WS, DP 555 BGRR, and BRS 369 RF obtained similar recommendations regarding adaptability parameters, while the others (TMG 41 WS, NUOPAL, DELTA OPAL, BRS 286, BRS 335, and BRS 368 RF) had divergent classifications between ECERSEM-AdaptStab and ER (Table 3). These differences can be attributed to the influence of environmental covariates, which are not considered by the traditional ER model. By capturing these environmental effects, ECERSEM-AdaptStab may generate different parameter estimates, potentially even leading to a difference in the classification of genotype suitability compared to ER.

Given these results, it can be seen that the new strategy brings with it the advantages of incorporating information from environmental covariates directly into the model and not just working with the averages of each location. Furthermore, it is possible to separate the impact of environmental covariates on both the genotypes and the index, with this impact being exactly associated with the experimental information (location and year) imposed on the model. However, like all regression models, the addition of parameters is totally related to the degrees of freedom of the model, which directly impacts a minimum of observations for the model to work.

Multiple studies, employing the same genotype data utilized in ECERSEM-AdaptStab, presented discrepancies between the ER approach and other model-based regression. These discrepancies involve both adaptability parameters and the coefficient of determination (R²). Nascimento et al. [25], using quantile and non-parametric regressions, revealed distinct adaptability patterns among genotypes TMG 41 WS, TMG 43 WS, BRS 286, and BRS 368 RF. Nascimento et al. [12] using a Bayesian segmented regression model found discrepancies with results obtained for NUOPAL and DELTA OPAL genotypes. Similarly, Teodoro et al. [31] also observed discordances with ER using the non-parametric method proposed by Lin and Binns [34] modified by Carneiro [35] for genotypes IMA CV 690, DP 555 BGRR, and BRS 369 RF. These contrasting classifications indicated the potential influence of methodological choices on adaptability and stability recommendations.

Selecting the most suitable method for analyzing adaptability and stability depends on the specific research objectives and the nature of the available data. For example, analysis of variance (ANOVA) methods allow to estimate the contributions of genotype, environment, and GEI to the overall performance variation [36,37]. Multi-environment trial analysis methods like AMMI and GGE Biplot combine analysis of variance (ANOVA) with other techniques like principal component analysis (PCA) to identify genotypes with both high mean yields and specific adaptation to certain environments [38]. Additionally, the regression-based methods utilize the analysis to assess stability and identify genotypes with high mean yields and low responsiveness to environmental changes [6,7]. Each approach has its own characteristics and can be used alone or as complementary to each other. However, unlike the ER method and other approaches, some alternative approaches do allow for the incorporation of environmental covariates (EC) into the process.

Similar to many other crops, EC significantly interacts with the genotype of cotton plants, affecting their overall production and potentially influencing their adaptability and stability. Among the EC, heat stress causes a potential decline in cotton yields [39,40,41]. Schlenker and Roberts [42] project that cotton-growing regions could face significant warming exceeding the global average. This warming trend has the potential to cause cotton yields to decline by up to 40% by 2100. In addition to heat stress, water stress caused significant reductions in days to first flower formation, plant height, number of bolls/plant, boll weight, and seed cotton yield [43,44,45]. Other studies have also shown that increased sunlight decreased the production of some cotton cultivars due to increased evaporation and water stress in some regions of China [46], but favored increased production of other cotton cultivars in other regions, an increase explained by favoring photosynthesis [47]. It was also found that wind speed in critical periods of the cotton production window can affect cotton sowing and seedling growth [48].

Considering our results, the genotypes TMG 41 WS, TMG 43 WS, IMA 5675 B2RF, IMA 08 WS, and BRS 368 RF were not affected directly by the EC in these specific season conditions (Figure 5 and Table S1). This result suggests that these genotypes are stable over the direct effects of EC variations. On the other hand, TS and T2M (related to the heat stress) affected six (IMA CV 690, NUOPAL, DP 555 BGRR, DELTA OPAL, BRS 286, and BRS 335) and five (IMA CV 690, NUOPAL, DP 555 BGRR, DELTA OPAL, and BRS 286) genotypes, respectively, being, in our study, the most important EC in this set of environments. It is important to note that if the number of seasons were different, the results would possibly change, as it could also be different in other models based on averages of the environments already proposed, as this would directly impact the value linked to the specific environment.

Specifically, TS negatively affected IMA CV 690 and DP 555 BGRR, while T2M and CLRSKY_SFC_PAR_TOT affected the same genotypes positively. IMA CV 690 was negatively affected by RH2M, while genotype DP 555 BGRR suffered a positive effect from PRECTOTCORR and WS2. While temperature (T2M) and light availability (CLRSKY_SFC_PAR_TOT) presented similar effects on these two genotypes, caution is warranted in their recommendation. Their differential responses to other EC, like humidity (RH2M), rainfall (PRECTOTCORR), and wind speed (WS2_MAX), highlight the importance of considering information about the EC. Therefore, the recommendation for genotype IMA CV 690 should be directed to regions that present milder temperatures and good CLRSKY_SFC_PAR_TOT conditions, but regions that do not present high RH2M should also be considered. On the other hand, genotype DP 555 BGRR should be recommended for regions that also have good CLRSKY_SFC_PAR_TOT conditions, but also with good precipitation conditions and without the occurrence of WS2_MAX. Genotype BRS 369 RF, contrary to genotype IMA CV 690, should be allocated to regions where the occurrence of higher RH2M is more frequent. Finally, BRS 368 RF was not significantly influenced by any of the EC used in this work; therefore, the recommendation would be to maintain the same as would traditionally be done in classical methods, as was observed by the ECERSEM-AdaptStab approach. These results showed that variations in EC have affected these genotypes, which potentially will change their recommendations depending on climatic variations.

The indirect effects of RH2M and PRECTOCORR on all genotypes through environmental index were not significant (Figure 6). Other environmental covariates showed consistently directional indirect effects. The low difference between these results is expected because they arise from the product of two parameters. The first one, gamma ($\mathsf{\gamma }$), represents the direct effect of each environmental covariate on the environmental index. The second parameter captures the direct effect of the environmental index on each genotype’s adaptability, which in this work ranges from 0.43 to 1.55. It is interesting to notice that some environmental covariates (ECs) did not present a direct effect on genotypes, but their influence was observed indirectly (Figure 6). This information, given by the SEM modeling proposal, is particularly valuable since even for ECs with no direct effect on genotypes, we can still uncover their influence by evaluating their indirect effect through the environmental index. Specifically, the genotypes TMG 41 WS, TMG 43 WS, IMA 5675 B2RF, IMA 08 WS, and BRS 368 RF presented indirect effects for TS, T2M, WS2M_MAX, and CLRSKY_SFC_PAR_TOT.

Regarding the total effects, it is possible to observe that TS and WS2M_MAX exhibited positive effects, while T2M, RH2M, PRECTOTCORR, and CLRSKY_SFC_TOT showed negative effects (Figure 7). This result allows the researcher to obtain responses about the total effect of the environmental covariate on the genotype. However, by partitioning this effect, it is possible to verify how these covariates directly influence genotypes or indirectly modulate their effects across the studied environments. In practice, by partitioning the total effects, the researchers improve their understanding of how the EC affects the genotypes, bringing more knowledge to help their decisions.

In theory of regression, including an additional explanatory variable (EC) can improve the model’s fit by capturing more variance in the trait [49]. This may also lead to smaller standard errors, making the coefficients of all variables more precise estimates. Additionally, omitting relevant variables can lead to several issues, including inflated model error and biased and inconsistent estimators, ultimately resulting in unreliable models [35]. However, introducing highly correlated explanatory variables can introduce multicollinearity. This can create unstable models with insignificant coefficients, making interpretation difficult [50]. The combined effects of multicollinearity and sample size limitations restrict the number of explanatory variables (ECs) we can effectively include in the model. The sample size is related to the choice of regression-based method for adaptability and stability analysis. To do that, the researcher should consider the number of environments included in the analysis. For instance, the Eberhart and Russell (ER) model requires at least three environments due to its two estimated parameters. Conversely, segmented regression models necessitate a minimum of six environments. These models typically allocate three environments below the average environmental index and three above to capture the non-linear response. As the proposed method incorporates EC, it necessitates fitting a model with more parameters. Researchers should account for this during the analysis.

Overall, ECERSEM-AdaptStab unveils the direct and indirect effects of environmental covariates on genotypes. Direct effects describe how individual ECs influence genotypes, while indirect effects capture how these factors combine to impact genotypes through an environmental index. In a practical setting, this information is very important for breeders as climate change intensifies. By understanding how genotypes respond directly and indirectly to environmental variations, breeders can make informed decisions about which lines to prioritize for breeding programs. Recent studies have explored the significant impact of climate change on crop breeding [51,52]. These studies emphasize the crucial role of crop breeding in adapting to these challenges, complementing crop management strategies and policy interventions to ensure global food security. However, to our knowledge, no existing adaptability and stability method explicitly incorporates information about ECs into its recommendation process.

5. Conclusions

In this study, a methodology based on SEM was proposed, which combines traditional parameters of adaptability and stability with the use of environmental covariates, proving to be advantageous given the new climate scenarios that are to come. Furthermore, a user-friendly shiny application (ECERSEM-AdaptStab) was created that incorporated all procedures, with the aim of facilitating use by plant breeders.

ECERSEM-AdaptStab identified that the genotypes TMG 41 WS, IMA CV 690, DP 555 BGRR, BRS 286 and BRS 369 RF should be recommended for favorable environments, while the genotypes TMG 43 WS, IMA 5675 B2RF, IMA 08 WS, NUOPAL, DELTA OPAL, BRS 335, and BRS 368 RF are more suitable for unfavorable environments.