Modeling temporal genetic variability using mixed models improves yield stability and selection efficiency in Coffea canephora

Abstract

Introduction:

Increasing climatic variability challenges Coffea canephora breeding programs to identify genotypes that combine high productivity with temporal stability across contrasting seasons.

Methods:

We evaluated 44 genotypes across four consecutive crop seasons (2022–2025) in eastern Minas Gerais, Brazil, using mixed linear models (REML/BLUP) with seven alternative variance-covariance structures for genetic and residual effects.

Results:

The flexible model M7 (unstructured genetic covariance matrix with year-specific residual variances) provided the best fit (lowest AIC and BIC). Plot-level heritability ranged from 0.64 to 0.68, genotype-mean heritability from 0.84 to 0.86, repeatability was 0.89, and selective accuracy ranged from 0.92 to 0.94. Genetic correlations revealed atypical behavior in 2023, driven by heat stress during grain filling. Relative selection efficiency increased cumulatively by 4.4% after four years. Genotypes Bicudo, A1, and AD1 combined the highest predicted genotypic values with elevated persistence indices.

Discussion:

Flexible mixed model approaches improve the reliability of genetic evaluation and support resilience-oriented selection strategies in C. canephora, enabling accelerated genetic gain and identification of superior genotypes adapted to variable cultivation conditions.

1 Introduction

Increasing climatic variability has intensified yield instability in coffee production systems, challenging the capacity of breeding programs to identify genotypes capable of sustaining high and stable productivity across contrasting years (Tavares et al., 2018; Ferreira et al., 2024). Across several producing regions, climate-driven heat and water stress have been associated with recurrent yield declines and marked interannual variability, with losses often reaching double-digit percentages (Pham et al., 2019). For instance, prolonged drought combined with recurrent temperature extremes in Brazil has reduced productive potential in key coffee-growing regions, leading to downward adjustments in 2025/26 output projections and contributing to price volatility in international markets¹. In perennial crops, where selection relies on repeated multi-year evaluations, such instability intensifies genotype × year interaction and increases the likelihood of recommending genotypes whose performance is contingent on transient environmental conditions (Fernandes Filho et al., 2025).

Coffea canephora, commonly known as Robusta or Conilon coffee, accounts for more than 40% of global coffee bean production². Its relevance stems from its tolerance to high temperatures, adaptation to low-altitude environments, and relative resistance to pests and diseases (Ferrão et al., 2024). Brazil is the second largest producer of the species and holds a substantial genetic heritage in germplasm collections, while continuously investing in breeding programs, thereby establishing itself as a global reference in the development of new cultivars (Jordaim et al., 2025; Adunola et al., 2024; Ferrão et al., 2020). Beyond direct physiological effects on plant metabolism and fruit development, climate change is also expected to alter the distribution and severity of pests and diseases, including the coffee leaf miner (Leucoptera coffeella) and fungal pathogens such as Hemileia vastatrix, further compounding yield instability and reinforcing the need for resilience-oriented selection strategies (Harelimana et al., 2022).

Brazilian breeding programs aim to sustainably increase productivity and production efficiency (Baqueta et al., 2024). However, in perennial crops such as C. canephora, increasing interannual climatic variability complicates the identification of genotypes that combine high productivity with resilience across contrasting seasons. Therefore, genotypes selected based on performance in a limited number of crop years may fail to sustain yield under subsequent climatic conditions, increasing the risk of productivity losses at the field level. Conversely, reducing the number of crop seasons to lower costs and accelerate selection cycles may compromise the accuracy of genetic estimates due to strong environmental effects (Adunola et al., 2023; Fernandes Filho et al., 2025). In this context, the estimation of genetic parameters such as repeatability and heritability is essential to balance cost, accuracy, and selection gain (Roka et al., 2024).

Traditionally, analysis of variance (i.e. ANOVA) and mixed linear models under REML/BLUP in their simplified forms have been employed to define the number of harvest years and to assess genotype-by-year interaction in coffee breeding (Mistro et al., 2008; Rocha et al., 2015). Nonetheless, these approaches assume homogeneous variances and restricted correlations across years, assumptions that are often unrealistic under variable climatic conditions. Such simplifications may obscure biologically meaningful differences in genotype performance and limit the capacity to correctly characterize temporal patterns of genetic stability. Longitudinal analyses of yield stability in C. canephora clones have confirmed that covariance structures such as CS, CSH, and UN provide adequate fit to multi-year data in this species (Cilas et al., 2011), and that temporal instability is a consistent feature of production even across extended evaluation periods of up to 14 years (Cilas et al., 2003). More recently, Covre et al. (2022) applied a Bayesian MCMC framework to 43 C. canephora genotypes evaluated across four harvests in two contrasting environments and reported broad-sense heritability of 0.28, highlighting the challenges of precise genetic estimation when temporal covariance structure is not explicitly modeled. In contrast, the use of mixed models incorporating alternative variance-covariance structures allows greater flexibility in modeling heterogeneous variances and non-uniform correlations enabling greater precision in the genetic evaluation of perennial crops such as coffee (Henderson, 1975; Patterson and Thompson, 1971; Piepho et al., 2008; Evangelista et al., 2023).

Despite these methodological advances, systematic comparisons among different variance-covariance structures in C. canephora remain scarce. Most analyses rely on simplified models, which restrict the understanding of the complex genotype-by-year interactions in this perennial species (Bergo et al., 2020). The adoption of more flexible structures, capable of accommodating heterogeneous variances and non-uniform correlations represents critical, yet still underexplored, step toward improving selection decisions in C. canephora breeding programs.

Therefore, understanding genotype-by-year interaction is crucial for identifying superior genotypes and determining the appropriate number of crop seasons required for reliable and resilient selection. Accordingly, the objectives of this study were to: (i) compare variance-covariance structures in C. canephora; (ii) estimate key genetic parameters relevant to yield stability and resilience; (iii) determine the minimum number of years required for reliable selection and resilient selection; and (iv) identify genotypes combining high productivity with persistence under variable climatic conditions in order to enhance breeding strategies for C. canephora in Brazil.

2 Materials and methods

2.1 Experimental area and cultivation conditions

The experimental area was located in the municipality of Aimorés, in the eastern region of Minas Gerais, Brazil, at 19°34′54.79″ S latitude, 41°23′00.99″ W longitude, and approximately 260 m altitude. The predominant soil type was classified as dystrophic Red-Yellow Latosol, according to the Brazilian Soil Classification System corresponding to the Oxisol in USDA Soil Taxonomy and Ferralsol in WRB (Santos et al., 2025).

The regional climate is tropical, characterized by hot and humid summers and dry winters, and classified as Aw under the Köppen system (Alvares et al., 2013; Beck et al., 2023). During the experimental period (2022-2025), the mean air temperature was 24.1 ± 3.9 °C (mean ± standard deviation), with absolute maximum and minimum temperatures of 36.1 °C and 13.0 °C, respectively. Total precipitation reached 862 mm, with an annual average of 216 mm. The mean relative humidity was 70.4 ± 16.8% (Figure 1)³. Notably, the 2023 growing season was characterized by recurrent high-temperature episodes during grain filling, providing a natural contrast for evaluating genotype performance under heat stress.

Figure 1

2.2 Genetic material and experimental design

A total of 44 Coffea canephora genotypes were evaluated, cultivated at a spacing of 3.2 m between rows and 0.8 m between plants (approximately 3,906 plants ha^-1), with two orthotropic stems per plant. The genotype set comprised three groups (Table 1): two Embrapa cultivars (category a), 23 promising Conilon genotypes originating from traditional producing regions in Espírito Santo (category c), and 19 seed derived selections from eastern Minas Gerais (category b named LMG). These LMG genotypes (category b) were initially identified as seed derived selections in farmers’ fields and were subsequently clonally propagated for inclusion in the trial. All genotypes were established simultaneously and evaluated at the same plant age. Detailed identification and classification of the evaluated genotypes is presented in Table 1.

The experiment was conducted using a randomized complete block design, with 44 treatments (genotypes) and three replications (blocks or experimental plots). Each plot consisted of five plants, with only the three central plants evaluated to minimize border effects. Crop management followed technical recommendations for coffee cultivation. Harvests were carried out annually from 2022 to 2025, totaling four consecutive crop seasons.

Productivity was determined based on the harvest of ripe fruits from each experimental plot, conducted separately by genotype. The harvested fruit volume was initially measured in liters per plot and subsequently converted to 60-kg bags of processed coffee per hectare. For this conversion, an equivalence of 320 liters of ripe fruits per 60-kg bag of processed beans was considered, with adjustments made according to plant density per hectare based on the adopted spacing.

2.3 Statistical analyses

Analyses were performed using the methodology of mixed linear models. The basic model considered was:

where represents the vector of phenotypic observations; is the overall mean associated with the vector of ones; is the vector of fixed effects of years with incidence matrix ; and is the vector of fixed effects of replications within years, associated with incidence matrix . The term corresponds to the vector of random genotypic effects, assumed as , where is the matrix of genetic variances and covariances across years, and is . The term is the vector of random permanent plot effects (), assumed as , where is the variance of permanent environmental effects and is the identity matrix of order 132 (number of plots), associated with incidence matrix . Finally, is the vector of random errors, assumed as , where is the matrix of residual variances across years and is the identity matrix of order (number of observations per year).

2.4 Modeling of genetic effects

The matrix was modeled with four alternative covariance structures, while two approaches were considered for the residual matrix , yielding seven models in total. These structures were defined a priori as a systematic progression from the most restrictive to the most flexible parameterisation applicable to repeated-measures data in plant breeding (Gilmour et al., 2009; Butler et al., 2023), covering the principal options available for modeling and . This progression was not derived from preliminary data screening but reflects the theoretical framework for longitudinal mixed models, in which progressively relaxing constraints on variances and correlations allows an empirical assessment of which assumptions are supported by the data.

Initially, a baseline model was fitted, adopting a homogeneous compound symmetry structure, which partitions the genetic variance into the main effect and the effect due to genotype-by-year interaction (GYI):

where is the number of genotypes (). This structure treats all crop seasons as exchangeable environments, assuming that the genetic expression of yield is equally variable and equally correlated across any pair of years, an assumption equivalent to a traditional repeatability model.

The second covariance structure was the diagonal (DIAG), which accounts for heterogeneity of genetic variances across years while assuming null paired covariances:

By allowing year-specific genetic variances while fixing all covariances at zero, DIAG acknowledges that the magnitude of genetic expression can differ across seasons but assumes that the relative ranking of genotypes is entirely independent between any two years, a biologically extreme assumption for a perennial crop evaluated under variable but related climatic conditions.

The next covariance structure was heterogeneous compound symmetry (CSH), which allows heterogeneous genetic variances across years while assuming a common correlation () between any pair of years:

CSH relaxes the variance constraint of CS while retaining a single correlation parameter shared across all year pairs, reflecting the assumption that the degree of genetic consistency between seasons is constant regardless of which two years are being compared.

Finally, to more flexibly account for heterogeneity of genetic variances, an unstructured (UN) model was fitted. This structure imposes no constraints on either variances or pairwise correlations, allowing the data to determine freely how genotypic expression varies across seasons and how consistently genotypes perform between any specific pair of years. It is the biologically most realistic option for perennial crops where each crop season integrates a unique combination of climatic conditions, bearing cycle stage, and cumulative plant developmental effects.

2.5 Modeling of residual effects

Two main approaches were considered for modeling residual effects. In the first case, homogeneity of residual variances across years was assumed, using the identity structure:

where is the common residual variance and is the identity matrix of order .

In a second step, models with heterogeneous residuals were fitted, allowing specific variances for each year. In this case, the residual matrix assumes a heterogeneous diagonal form:

where represents the residual variance estimated for the j-th year.

2.6 Model selection

Model fit was evaluated using the Akaike Information Criterion (AIC; Akaike, 1974) and the Bayesian Information Criterion (BIC; Schwarz, 1978). The AIC is defined as

where denotes the maximum value of the (restricted) likelihood function and is the number of estimated parameters. The BIC was computed as

with representing the number of observations. In both cases, smaller values indicate a more parsimonious balance between model fit and complexity.

2.7 Estimation of genetic and non-genetic parameters

The estimation of variance components and genetic parameters was performed using mixed models, considering both the baseline model (M1) and the best-fitting model. From the variance components obtained, several genetic parameters of interest were calculated.

The individual phenotypic variance was defined as the sum of genetic variance, genotype-by-year interaction variance, permanent environmental variance, and residual variance:

The repeatability coefficient was estimated as the ratio between the sum of genetic and permanent environmental variances and the total phenotypic variance:

Plot-level heritability was calculated as the ratio between genetic variance and phenotypic variance:

Genotype-mean heritability was obtained by considering the number of years () and the number of replications (), according to the following expression:

Cullis heritability (Cullis et al., 2006) was estimated from the average variance of prediction errors () of pairwise genotype differences:

Selective accuracy was obtained following Mrode (2014), calculated as:

where PEV represents the prediction error variance associated with each genotype. Parameters were estimated both in the baseline model and in the best-fitting model for comparative purposes. When the best-fitting model presented a heterogeneous structure, values of and were defined separately for each year.

Genetic correlations across years were derived from the variance-covariance matrices of genotypic effects fitted under different structures. In the case of the unstructured (UN) model, the genetic correlation between years and was estimated as the ratio between the genetic covariance and the product of the square roots of the marginal genetic variances:

where represents the genetic covariance between years and , and and are the specific genetic variances for each year.

Relative selection efficiency was estimated based on the cumulative repeatability of yield across years, following the approach proposed by Resende (2002). For this purpose, the genotypic variance-covariance matrix across years was constructed, from which the average repeatability was obtained. Efficiency was then calculated as:

where represents the number of years considered and the cumulative repeatability.

Yield persistence of genotypes was estimated from the predicted genotypic values () obtained from the best-fitting model (Rocha et al., 2018). For each year (), the ideotype was defined as the highest BLUP observed (). Subsequently, for each genotype (), the sum of squared differences between its predicted value and the ideotype across years was calculated. The persistence measure was obtained as the ratio between this sum and the total accumulated across all genotypes, expressed on a relative scale:

where corresponds to the number of years evaluated and to the number of genotypes. Higher values of indicate greater proximity of the genotype to the ideotype across years, reflecting higher yield stability and persistence.

All analyses were implemented in the RStudio integrated development environment (IDE), using the R programming language (version 4.5.1; R Core Team, 2025). The ASReml-R package (Butler et al., 2023), version 4.2.0.355, was employed for model fitting, while ggplot2 (Wickham, 2016) was used for graphical visualization. A fully reproducible R script detailing data processing, mixed model specification, model updating, and the computation of all derived parameters is provided in the Supplementary Material.

3 Results

3.1 Model selection

Model comparison revealed clear differences in goodness of fit among the seven variance–covariance structures evaluated (Table 2). The model combining an unstructured genetic covariance matrix with heterogeneous residual variances (M7) provided the best fit according to both AIC and BIC criteria and was therefore selected for subsequent analyses (Table 2). Simpler models assuming homogeneous genetic variances and correlations performed less well. Model M1, which employed compound symmetry for genetic effects and homogeneous independent residuals, yielded an AIC of 4190.61 and a BIC of 4258.92, with only three parameters estimated, serving as a reference equivalent to a traditional ANOVA model. Allowing heterogeneity of genetic variances across years (M2 and M3) did not substantially improve model fit, although M3 showed a slight improvement relative to M2.

Models allowing greater flexibility in genetic covariance structure consistently outperformed simpler alternatives. The unstructured genetic model with homogeneous residuals (M4) substantially reduced AIC values (4176.91), while the inclusion of heterogeneous residual variances (M7) resulted in the lowest AIC among all models tested (4173.71). This difference relative to M4 indicates that modeling residual heterogeneity contributed to a more adequate fit. Among models with heterogeneous compound symmetry for , M6 provided a better fit than M3, but still inferior to the unstructured models. Variance components and genetic parameters estimated under all seven models are provided in Supplementary Table 1.

Despite variations in information criteria, predictive accuracy was consistently high across all models, ranging from 0.91 to 0.93. This indicates that model choice primarily affected variance partitioning and biological interpretation rather than the stability of genotype rankings. Considering the balance between model fit and flexibility, M7 was selected as the most appropriate model for subsequent analyses, and all genetic parameters reported hereafter refer to estimates obtained under this model.

3.2 Genetic and non-genetic parameters

Estimates of variance components and genetic parameters obtained from the baseline model (M1) and the best-fitting model (M7) are presented in Table 3. The simplified covariance structure assumed in M1 resulted in relatively low genetic variance () and a high contribution of genotype-by-year interaction (), reflecting the limitation of this model in adequately separating genetic effects from temporal variation. In contrast, the more flexible model M7, which allowed year-specific genetic variances and heterogeneous residuals, revealed substantially higher estimates, ranging from 822.20 in 2024 to 1319.03 in 2022. Residual variances were also heterogeneous across years. This improved variance partitioning resulted in higher phenotypic variance estimates and greater consistency of genetic parameters across crop seasons.

Plot-level heritability () was low under the baseline M1 (0.19), but increased markedly under M7, reaching values between 0.64 and 0.68 across years. Similarly, genotype-mean heritability () was also higher in the best-fitting model (0.84–0.86), compared with 0.55 in the baseline model. Cullis heritability () remained high under both models but showed greater stability across years under M7 (0.85–0.89). The repeatability coefficient () increased from 0.56 in M1 to 0.89 ( cumulative) in M7, highlighting a much higher consistency of genotype performance across years when temporal heterogeneity was explicitly modeled.

Repeatability showed consistently high values across the four years of evaluation, ranging from 0.84 to 0.87, thereby indicating stable expression of yield across crop cycles (Figure 2A). Relative selection efficiency increased steadily with the addition of evaluation years, rising from 1.00 in the first year (2022) to 1.04 in the fourth year (2025), representing a cumulative gain of approximately 4.4% (Figure 2B). Selective accuracy remained high across all years, oscillating between 0.92 and 0.94, and confirming the reliability of the predicted genotypic values (Figure 2C).

Figure 2

The analysis of genetic correlations among harvest years revealed wide variation in the magnitude of coefficients (Figure 3). The strongest associations were observed between 2022 and 2024 (0.70) and between 2024 and 2025 (0.64), indicating greater genetic stability in these pairs of environments. In contrast, the correlation between 2022 and 2025 was low (0.21), suggesting limited consistency and reordering of genotypes. The 2023 season, characterized by extreme temperatures during grain filling, exhibited markedly reduced genetic correlations with all other years (ranging from -0.10 to 0.21), characterizing atypical behavior and strong genotype-by-environment interaction during that season.

Figure 3

These results indicate that 2022, 2024, and 2025 share greater genetic similarity, whereas 2023 diverges substantially from the observed pattern. From a biological standpoint, this divergence reflects differential genotypic sensitivity to the climatic conditions of that season: recurrent heat episodes above 35 °C during grain filling in 2023 likely imposed asymmetric stress across genotypes, amplifying differences in physiological buffering capacity and producing a reordering of performance ranks that no model assuming stable inter-year correlations could detect (Hasanuzzaman et al., 2013; de Oliveira et al., 2020). The detection of this anomalous season is itself a direct product of the unstructured covariance approach, under any model constraining correlations to be equal across year pairs, the atypical behavior of 2023 would have been absorbed into a pooled GxY term, obscuring rather than characterizing the interaction.

3.3 Selection efficiency and yield persistence

The joint analysis of predicted genotypic values for yield and the persistence index enabled the identification of superior genotypes in C. canephora (Figure 4). Genotypes Bicudo, A1 and AD1, which belong to the group of promising Conilon genotypes traditionally cultivated in Espírito Santo, stood out by simultaneously exhibiting the highest predicted genotypic values (121, 119 and 118 bags ha^-1, respectively) and elevated persistence indices, demonstrating consistent performance across years. Other genotypes, such as LMG1 and LMG3, which represent new seed-derived selections from eastern Minas Gerais, as well as Valcir P and LB1, both recognized as promising genotypes from Espírito Santo, also ranked among the most promising by combining satisfactory yield with stable performance. In contrast, genotypes LMG10, LMG14 and LMG7, all originating from the group of new seed-derived selections, showed low values for both criteria, reflecting inferior performance and reduced production regularity. These findings underscore the usefulness of jointly evaluating predicted genotypic values and persistence indices for selecting genotypes that integrate high productivity with stable performance across contrasting crop seasons.

Figure 4

4 Discussion

Repeatability studies are fundamental in the breeding of perennial species such as coffee as selection decisions rely on repeated measurements of the same genotypes over time. These long-term evaluations inevitably expose plant material to fluctuating environmental conditions, leading to heterogeneous variances and complex covariance patterns driven by genotype × environment interactions operating at multiple biological scales (Chaves et al., 2023). Under such conditions, the use of statistical models capable of accommodating this complexity becomes essential to ensure biologically meaningful interpretation and reliable genetic parameter estimation.

The comparison among variance-covariance structures underscores the need for more flexible models to capture year-to-year dynamics in C. canephora. The inferior performance of M2 and M3 models (Table 2), even with the inclusion of genetic heterogeneity, suggests that simply modeling distinct variances is insufficient to adequately represent temporal correlation. In contrast, unstructured models yielded improved statistical fit, demonstrating that genotype-by-environment interaction in perennial species rarely follows simplified patterns of homogeneity (Alves et al., 2023).

Among the models evaluated, M4 and M7 stood out, both incorporating an unstructured . The additional improvement in fit observed in M7, which combined an unstructured with a diagonal , shows that accounting for heterogeneous residuals contributes to a more accurate description of experimental variability. Models with heterogeneous compound symmetry for , such as M6, performed at an intermediate level, better than M3 but still inferior to unstructured models. Although predictive accuracy remained high across all models, the choice of M7 is justified by its superior statistical fit, greater biological consistency, and robustness for downstream inference. Importantly, this indicates that model choice primarily affects variance partitioning and biological interpretation rather than altering genotype rankings per se.

The contrast between models reveals the influence of covariance structure on the estimation of genetic parameters (Table 3). Model M1, by assuming homogeneity across years, substantially underestimated genetic variance, resulting in very low plot-level heritability () and only moderate genotype-mean heritability (). This behavior reflects a dilution of the genetic signal when temporal environmental variation is inadequately modeled. In contrast, model M7, by allowing year-specific genetic variances and heterogeneous residuals, substantially increased plot-level heritability (0.64–0.68) and genotype-mean heritability (0.84–0.86), while maintaining consistently high (0.85–0.89). This behavior indicates that more flexible modeling captures the data structure more realistically, reflecting the complexity of genotype-by-environment interaction.

The biological expectation of heterogeneous variances and unstructured correlations in perennial crop systems follows directly from the nature of multi-year field evaluations. Each crop season integrates a unique combination of rainfall distribution, temperature regime, bearing cycle stage, and cumulative plant developmental history, all of which affect fruit set, grain expansion, and carbon partitioning differently across genotypes. Under these conditions, assuming that genetic variance is constant across years, as in M1 and M2, or that any two seasons are equally correlated, as in CS and CSH, is biologically indefensible, not merely statistically suboptimal. The variation in genetic variance observed across years under M7 (σ²g ranging from 822.2 in 2024 to 1319.0 in 2022) reflects precisely this: in seasons where environmental conditions amplify differences in physiological buffering capacity among genotypes, the genetic signal is stronger and more reliably separated from residual noise. Longitudinal analyses in C. canephora have confirmed that flexible covariance structures provide better fit to multi-year yield data than restricted models (Cilas et al., 2011), and that the magnitude of genetic variance itself varies meaningfully across seasons in this species (Cilas et al., 2003).

The coherence between heritability parameters and genetic correlations reinforces this interpretation (Figure 3). Years 2022, 2024, and 2025 exhibited moderate to high correlations, reflecting greater consistency in genotype performance and justifying the elevated heritability estimates observed in M7. In contrast, 2023 showed very low correlations with the other years, signaling strong environmental influence and highlighting the need for models that account for heterogeneous residuals. The near-perfect repeatability in M7 () and consistently high selective accuracy (), as shown in Figure 2, confirm that the prediction of genetic values is robust when the variance structure is properly specified.

The consistently high repeatability values observed across the four years of evaluation indicate that a large proportion of phenotypic variation is explained by permanent genetic effects, conferring greater consistency to genotype performance across crop seasons (Figure 3). According to Resende and Alves (2022), repeatability coefficients above 0.80 enable reducing the number of harvests required for reliable selection without compromising the accuracy of estimates. This result is corroborated by the progressive increase in relative selection efficiency, which reached a cumulative gain of 4.4% in the fourth year. Such behavior underscores the importance of multiple crop seasons to consolidate the prediction of genetic values, as also reported in studies of perennial species (Chaves et al., 2022; Dessauw et al., 2024).

Selective accuracy remained high in all years (0.92–0.94; Figure 2C), confirming the robustness of the estimates obtained and the reliability of the selection process. Accuracy values above 0.90 are considered highly precise and ensure that the selection reflects the true genetic merit of individuals (Cullis et al., 2006; Resende, 2016). These results demonstrate that, when the variance structure is properly specified, it is possible to shorten the evaluation period without loss of reliability, thereby increasing the efficiency of the breeding program.

Moreover, the robustness of the estimates obtained in this study derives not only from the more flexible statistical modeling (M7) but also from the quality of the balanced experimental design adopted. The use of randomized complete blocks contributed to reducing uncontrolled environmental variability, ensuring fairer comparisons among genotypes and decreasing residual variance. This experimental structure favored the attainment of more consistent heritability estimates, high repeatability, and superior selective accuracy, as observed in the results. According to Piepho et al. (2008) and Resende (2016), balanced designs maximize statistical efficiency, enhance the precision of genetic estimates, and strengthen the reliability of predictions. While longer evaluation periods may further refine parameter estimates, the present results demonstrate that robust and reliable selection can already be achieved within four crop seasons when appropriate modeling is applied.

The pattern of genetic correlations across years revealed marked heterogeneity in genotype performance, with 2023 behaving as an atypical season (Figure 3). The near-zero or negative correlations between 2023 and other years indicate strong genotype × environment interaction, while moderate to high correlations among 2022, 2024, and 2025 suggest greater genetic consistency during these periods. This atypical behavior is consistent with the climatic conditions recorded in 2023, when recurrent episodes of heat stress occurred, with temperatures exceeding 35 °C during the critical period of grain expansion and filling. Under such conditions, the photosynthetic apparatus may undergoes accelerated degradation, with damage to thylakoid protein complexes and reduced efficiency of photosystem II (Berry and Bjorkman, 1980; Hüve et al., 2011). Simultaneously, elevated temperatures induce stomatal closure as a protective mechanism against excessive water loss, thereby compromising CO₂ assimilation and limiting the production of photoassimilates precisely during the phase of highest metabolic demand by the fruits (Hasanuzzaman et al., 2013). Similarly, in coffee the elevated temperatures impact metabolism, fruit production and quality (Läderach et al., 2017; de Oliveira et al., 2020; Thioune et al., 2020; Ahmed et al., 2021), which could explain the results in 2023. A Bayesian MCMC analysis of 43 C. canephora genotypes evaluated across four harvests in two contrasting environments reported broad-sense heritability of 0.28 (Covre et al., 2022), a value substantially lower than those obtained in the present study under M7 (0.64–0.68). This contrast suggests that the explicit modeling of temporal covariance structure within the REML/BLUP framework yields more precise partitioning of genetic variance. Covre et al. (2022) also noted that climatic variables would be strong candidates as covariates in future models, an observation that aligns with the atypical behavior of 2023 documented here and points toward a productive direction for subsequent analyses.

The joint analysis of mean predicted genotypic values and yield persistence enabled the identification of genotypes that combine high performance with greater stability across years (Figure 4). Genotypes such as Bicudo, A1 and AD1, all of which belong to the group of promising Conilon genotypes traditionally cultivated in Espírito Santo, consistently ranked among the most favorable, aligning with the high heritability and repeatability estimates obtained under M7. Persistence therefore captures a resilience-related property, reflecting the capacity of genotypes to maintain productive performance despite strong interannual environmental fluctuations.

This approach becomes particularly relevant when confronted with the genetic correlation patterns observed over years, especially the atypical behavior of 2023, marked by strong genotype × environment interaction. High-persistence genotypes (Bicudo, A1 and AD1) exhibited reduced sensitivity to extreme conditions, maintaining relatively stable rankings across crop cycles. The ability of a genotype to maintain its relative rank even under heat stress, as recorded in 2023, reflects physiological resilience and greater predictability of productive behavior, attributes essential for ensuring consistent economic returns to producers (Malosetti et al., 2013). This stability contrasts with materials that, although showing good performance under specific conditions, experienced drastic shifts in classification across crop cycles. In contrast, genotypes such as LMG10, LMG14 and LMG7, which represent new seed-derived selections identified in eastern Minas Gerais, in addition to presenting low predicted genotypic means, also exhibited low persistence and greater variability in correlations across years. These genotypes tend to show erratic performance, hindering reliable recommendations and increasing the risk of failure in commercial production systems.

From a practical standpoint, the selection of genotypes that integrate high-predicted yield with elevated persistence helps reduce risks associated with genotype × environment interaction and enhances breeding efficiency (Cullis et al., 2006; Resende and Alves, 2022). The integration of genetic parameters, year-to-year correlations and persistence indices provides a solid basis for safer and more sustainable selection decisions, prioritizing not only absolute productivity but also the consistency of productivity across multiple cycles and under diverse climatic scenarios. This strategy aligns with the growing need to develop resilient cultivars capable of maintaining stable yield in the face of increasingly frequent and intense climate variability.

Some limitations of this study merit explicit acknowledgement. The experiment was conducted at a single location in eastern Minas Gerais, which limits the generalizability of the genotype rankings to other production environments. The four-year evaluation window, while sufficient to demonstrate the advantages of flexible modeling and to detect the atypical behavior of 2023, does not allow definitive conclusions about long-term yield stability, longer series would be required to determine whether the patterns of genetic correlation observed here are consistent over extended periods, as evaluated for up to 14 years in diallel populations of C. canephora by Cilas et al. (2003). Multi-environment validation would further clarify whether genotypes such as Bicudo, A1, and AD1, identified as high-yield and high-persistence at Aimorés, maintain these attributes across the broader range of edaphoclimatic conditions encountered in Brazilian C. canephora production regions.

5 Conclusions

Flexible variance-covariance models, particularly the unstructured model with heterogeneous residuals (M7), yielded superior statistical fit and a more realistic representation of genotype-by-year interaction in C. canephora. By allowing year-specific genetic variances and covariances, this model produced more accurate estimates of heritability, repeatability, and selective accuracy, demonstrating that reliable selection can be conducted within shorter evaluation cycles when appropriate modeling is applied.

The atypical genetic correlations observed in 2023, driven by heat stress during grain filling, underscore the need for multiple crop seasons to reduce the risk of biased selection decisions, a conclusion that would have been inaccessible under models assuming homogeneous covariance structures. The joint analysis of predicted genotypic values and the persistence index identified Bicudo, A1, and AD1 as genotypes combining high productivity with stable performance across contrasting seasons.

These results are subject to the limitation of a single-location design and a four-year evaluation window. Multi-environment and longer-term validation would strengthen the generalizability of the genotype rankings identified here and clarify whether the climatic drivers of the 2023 anomaly operate consistently across the broader range of C. canephora production environments in Brazil. Incorporating climatic covariates into the modeling framework represents a natural direction for future work.

References

• 1

  AdunolaP.FerrãoM. A. G.FerrãoR. G.FonsecaA. F. A.VolpiP. S.ComérioM.et al. (2023). Genomic selection for genotype performance and environmental stability in Coffea canephora. G3 (Bethesda)13, 6. doi: 10.1093/g3journal/jkad062. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AdunolaP.FloresE. T.Riva-SouzaE. M.FerrãoM. A. G.SenraJ. F. B.ComérioM.et al. (2024). A comparison of genomic and phenomic selection methods for yield prediction in Coffea canephora. Plant Phenome J.7, 1. doi: 10.1002/ppj2.20109. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  AhmedS.BrinkleyS.SmithE.SelaA.TheisenM.ThibodeauC.et al. (2021). Climate change and coffee quality: Systematic review on the effects of environmental and management variation on secondary metabolites and sensory attributes of Coffea arabica and Coffea canephora. Front. Plant Sci.12, 708013. doi: 10.3389/fpls.2021.708013. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  AkaikeH. (1974). A new look at the statistical model identification. IEEE Trans. Automat. Contr.19, 716–723. doi: 10.1109/TAC.1974.1100705. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  AlvaresC. A.StapeJ. L.SentelhasP. C.GonçalvesJ. L. M.SparovekG. (2013). Köppen’s climate classification map for Brazil. Meteorol. Z.22, 711–728. doi: 10.1127/0941-2948/2013/0507. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  AlvesA. K. S.ChavesS. F. S.AraújoM. S.MalikouskiR. G.AlmeidaC. M. V. C.DiasL. A. S. (2023). Improving multi-harvest data analysis in cacao breeding using random regression. Euphytica220, 1. doi: 10.1007/s10681-023-03270-6. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  BaquetaM. R.DinizP. H. G. D.PereiraL. L.AlmeidaF. L. C.ValderramaP.PalloneJ. A. L. (2024). An overview on the Brazilian Coffea canephora scenario and the current chemometrics-based spectroscopic research. Food Res. Int.194, 114866. doi: 10.1016/j.foodres.2024.114866. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  BeckH. E.McVicarT. R.VergopolanN.BergA.LutskoN. J.DufourA.et al. (2023). High-resolution (1 km) Köppen–Geiger maps for 1901–2099 based on constrained CMIP6 projections. Sci. Data10, 724. doi: 10.1038/s41597-023-02549-6. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  BergoC. L.MiqueloniD. P.LunzA. M. P.AssisG. M. L. (2020). Estimation of genetic parameters and selection of Coffea canephora progenies evaluated in Brazilian Western Amazon. Coffee Sci.15, 1–10.

  • Google Scholar

• 10

  BerryJ.BjorkmanO. (1980). Photosynthetic response and adaptation to temperature in higher plants. Annu. Rev. Plant Physiol.31, 491–543. doi: 10.1146/annurev.pp.31.060180.002423. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  ButlerD. G.CullisB. R.GilmourA. R.GogelB. G.ThompsonR. (2023). ASReml-R reference manual Vol. 4.2 (Hemel Hempstead, UK: VSN International Ltd).

  • Google Scholar
  • View reference in article

• 12

  ChavesS. F. S.AlvesR. S.DiasL. A. S.AlvesR. M.DiasK. O. G.EvangelistaJ. S. P. C. (2023). Analysis of repeated measures data through mixed models: An application in Theobroma grandiflorum breeding. Crop Sci.63, 2131–2144. doi: 10.1002/csc2.20995. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  ChavesS. F. S.DiasL. A. S.AlvesR. S.AlvesR. M.JoséA. R. M.AlmeidaC. M. V. C. (2022). Number of harvest years and selection for productivity, witches’ broom resistance, stability, and adaptability in cacao. Agron. J.114, 3234–3245. doi: 10.1002/agj2.21149. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  CilasC.BouharmontP.Bar-HenA. (2003). Yield stability in Coffea canephora from diallel mating designs monitored for 14 years. Heredity91, 528–532. doi: 10.1038/sj.hdy.6800351. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  CilasC.MontagnonC.Bar-HenA. (2011). Yield stability in clones of Coffea canephora in the short and medium term: longitudinal data analyses and measures of stability over time. Tree Genet. Genomes7, 421–429. doi: 10.1007/s11295-010-0344-4. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  CovreA. M.da SilvaF. A.OliosiG.CorreaC. C. G.VianaA. P.PartelliF. L. (2022). Multi-environment and multi-year Bayesian analysis approach in Coffee canephora. Plants11, 3274. doi: 10.3390/plants11233274. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  CullisB. R.SmithA. B.CoombesN. E. (2006). On the design of early generation variety trials with correlated data. J. Agric. Biol. Environ. Stat.11, 381–393. doi: 10.1198/108571106X154443

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  de OliveiraR. R.RibeiroT. H. C.CardonC. H.FedeniaL.MaiaV. A.BarbosaB. C. F.et al. (2020). Elevated temperatures impose transcriptional constraints and elicit intraspecific differences between coffee genotypes. Front. Plant Sci.11, 1113. doi: 10.3389/fpls.2020.01113. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  DessauwD.Phillips-MoraW.Mata-QuirósA.BastideP.JohnsonV.Castillo-FernándezJ.et al. (2024). Temporal behaviour of cacao clone production over 18 years. Agron. Sustain. Dev.44, 3. doi: 10.1007/s13593-024-00967-3. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  EvangelistaJ. S. P. C.PeixotoM. A.CoelhoI. F.FerreiraF. M.de Souza MarçalT.AlvesR. S.et al. (2023). Modeling covariance structures and optimizing Jatropha curcas breeding. Tree Genet. Genomes19, 21. doi: 10.1007/s11295-023-01596-9. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  Fernandes FilhoC. C.Lima BarriosS. C.SantosM. F.NunesJ. A. R.do ValleC. B.JankL.et al. (2025). Assessing genotype adaptability and stability in perennial forage breeding trials using random regression models for longitudinal dry matter yield data. G3 (Bethesda)15, 3. doi: 10.1093/g3journal/jkae306. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  FerrãoR. G.FerrãoM. A. G.VolpiP. S.FonsecaA. F. A.Verdin FilhoA. C.ComérioM. (2020). Cultivares de café Conilon e Robusta. Inf. Agropecu.41, 17–25.

  • Google Scholar
  • View reference in article

• 23

  FerrãoM. A. G.Riva-SouzaE. M.AzevedoC.VolpiP. S.FonsecaA. F. A.FerrãoR. G.et al. (2024). Robust and smart: Inference on phenotypic plasticity of Coffea canephora reveals adaptation to alternative environments. Crop Sci. doi: 10.1002/csc2.21298. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  FerreiraM. L.von dos Santos VelosoR.de OliveiraG. S.QueirozR. B.AraújoF. H. V.de AndradeA. M.et al. (2024). Effects of the climate change scenario on Coffea canephora production in Brazil using modeling tools. Trop. Ecol.65, 559–571. doi: 10.1007/s42965-024-00350-z. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  GilmourA. R.GogelB. G.CullisB. R.WelhamS. J.ThompsonR. (2009). ASReml user guide release 3.0 (Hemel Hempstead, UK: VSN International Ltd).

  • Google Scholar
  • View reference in article

• 26

  HarelimanaA.RukazambugaD.HanceT. (2022). Pests and diseases regulation in coffee agroecosystems by management systems and resistance in changing climate conditions: a review. J. Plant Dis. Prot.129, 1041–1052. doi: 10.1007/s41348-022-00628-1. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  HasanuzzamanM.NaharK.AlamM.RoychowdhuryR.FujitaM. (2013). Physiological, biochemical, and molecular mechanisms of heat stress tolerance in plants. Int. J. Mol. Sci.14, 9643–9684. doi: 10.3390/ijms14059643. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  HendersonC. R. (1975). Best linear unbiased estimation and prediction under a selection model. Biometrics. 31, 423–447. doi: 10.2307/2529430

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  HüveK.BicheleI.RasulovB.NiinemetsU. (2011). When it is too hot for photosynthesis: heat-induced instability of photosynthesis in relation to respiratory burst, cell permeability changes and H₂O₂ formation. Plant Cell Environ.34, 113–126. doi: 10.1111/j.1365-3040.2010.02229.x. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  JordaimR. B.ColodettiT. V.RodriguesW. N.SallesR. A.AmaralJ. F. T.MacielL. S.et al. (2025). Genotypic performance of Coffea canephora at transitional altitudes for climate-resilient coffee cultivation. Horticulturae11, 595. doi: 10.3390/horticulturae11060595. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  LäderachP.Ramirez–VillegasJ.Navarro-RacinesC.ZelayaC.Martinez–ValleA.JarvisA. (2017). Climate change adaptation of coffee production in space and time. Clim. Change141, 47–62. doi: 10.1007/s10584-016-1788-9. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  MalosettiM.RibautJ.-M.van EeuwijkF. A. (2013). The statistical analysis of multi-environment data: modeling genotype-by-environment interaction and its genetic basis. Front. Physiol.4, 44. doi: 10.3389/fphys.2013.00044. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  MistroJ. C.FazuoliL. C.Guerreiro FilhoO.SilvarollaM. B.Toma-BraghiniM. (2008). Determination of the number of years in Arabic coffee progenies selection through repeatability. Crop Breed. Appl. Biotechnol.8, 79–84. doi: 10.12702/1984-7033.v08n01a11

  • CrossRef
  • Google Scholar

• 34

  MrodeR. A. (2014). Linear models for the prediction of animal breeding values (Wallingford, UK: CABI).

  • Google Scholar

• 35

  PartelliF. L.LouzadaL. P.OliosiG.CampanharoA.CovreA. M.AlbertoN. J.et al. (2024). Research and development in Conilon and Robusta coffee (São Mateus, ES: Khas Editora).

  • Google Scholar

• 36

  PattersonH. D.ThompsonR. (1971). Recovery of inter-block information when block sizes are unequal. Biometrika58, 545–554. doi: 10.2307/2334389

  • CrossRef
  • Google Scholar

• 37

  PhamY.Reardon-SmithK.MushtaqS.DeoR.Nguyen-HuyT.StoneR.et al. (2019). The impact of climate change and variability on coffee production: a systematic review. Clim. Change156, 609–630. doi: 10.1007/s10584-019-02538-y. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  PiephoH.-P.MöhringJ.MelchingerA. E.BüchseA. (2008). BLUP for phenotypic selection in plant breeding and variety testing. Euphytica161, 209–228. doi: 10.1007/s10681-007-9449-8. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  R Core Team (2025). R: A language and environment for statistical computing (Vienna: R Foundation for Statistical Computing).

  • Google Scholar

• 40

  ResendeM. D. V. (2002). Genética biométrica e estatística no melhoramento de plantas perenes (Brasília: Embrapa Informação Tecnológica).

  • Google Scholar

• 41

  ResendeM. D. V. (2016). Software Selegen-REML/BLUP: a useful tool for plant breeding. Crop Breed. Appl. Biotechnol.16, 330–339. doi: 10.1590/1984-70332016v16n4a49. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  ResendeM. D. V.AlvesR. S. (2022). Statistical significance, selection accuracy, and experimental precision in plant breeding. Crop Breed. Appl. Biotechnol.22, e42712238. doi: 10.1590/1984-70332022v22n3a31. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  RochaJ. R. A. S. C.MarçalT. S.SalvadorF. V.SilvaA. C.MaChadoJ. C.CarneiroP. C. S. (2018). Genetic insights into elephantgrass persistence for bioenergy purpose. PloS One13, e0203818. doi: 10.1371/journal.pone.0203818. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  RochaR. B.RamalhoA. R.TeixeiraA. L.SouzaF. F.CruzC. D. (2015). Adaptability and stability of Coffea canephora coffee bean yield. Cienc. Rural45, 1531–1537. doi: 10.1590/0103-8478cr20141554. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  RokaP.ShresthaS.AdhikariS. P.NeupaneA.ShreepailiB.BistaM. K. (2024). A review on genetic parameters estimation, trait association, and multivariate analysis for crop improvement. Arch. Agric. Environ. Sci.9, 618–625. doi: 10.26832/24566632.2024.0903029

  • CrossRef
  • Google Scholar

• 46

  SantosH. G.JacomineP. K. T.AnjosL. H. C.OliveiraV. A.LumbrerasJ. F.CoelhoM. R.et al. (2025). Sistema Brasileiro de Classificação de Solos (Brasília: Embrapa).

  • Google Scholar

• 47

  SchwarzG. (1978). Estimating the dimension of a model. Ann. Stat.6, 461–464. doi: 10.1214/aos/1176344136

  • CrossRef
  • Google Scholar

• 48

  TavaresP. S.GiarollaA.ChouS. C.SilvaA. P.LyraA. A. (2018). Climate change impact on the potential yield of arabica coffee in southeast Brazil. Reg. Environ. Change18, 873–883. doi: 10.1007/s10113-017-1236-z. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 49

  ThiouneE.StricklerS.GallagherT.DuitamaJ.TorresJ. C.QuinteroC.et al. (2020). Temperature impacts the response of Coffea canephora to decreasing soil water availability. Trop. Plant Biol.13, 236–250. doi: 10.1007/s12042-020-09254-3. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 50

  WickhamH. (2016). ggplot2: Elegant graphics for data analysis (New York: Springer).

  • Google Scholar