The analyzed datasets correspond to three maize trials established at the Brazilian Agricultural Research Corporation (EMBRAPA) Maize and Sorghum headquarters in Sete Lagoas city at the Minas Gerais state, Brazil (19°28′S, 44°15′W). These trials were conducted in three consecutive agronomical years (2014/2015, 2015/2016, and 2016/2017), and each agronomical year was considered as a different environment (E1 = 2014/2015 year; E2 = 2015/2016 year, and E3 = 2016/2017 year). Genotypes were planted in a lattice design, with two repetitions. In total 373 single-cross hybrids were evaluated: 146 in environment 1 (E1), 145 in E2, and 150 in E3. Furthermore, care was taken to allow connectivity between environments: 33 hybrids are common between E1 and E2, 15 between E1 and E3, 37 between E2 and E3, and 13 across the three environments (Figure 1).

All 373 single-cross hybrids were scored for fumonisin concentration in parts per million (ppm) (FUMO), and this information was missing for the parents of these hybrids. The determination of the FUMO trait was initially based on a 500 g sample of corn grains that required a meticulous quantification process. The quantification took place at the Laboratory of Food Safety at Embrapa Milho e Sorgo where a 10 g subsample was finely ground. Then fumonisins were extracted in a solution of 100 mL water/methanol mixture (20/80), and 5 g NaCl using a blender for 1 min. The resulting mixture was then filtered through Whatman paper, and an aliquot of 10 mL of the filtered extract was subsequently diluted in 40 mL of 0.1% phosphate Tween-20 solution (phosphate buffer). This solution was then filtered again using a 1.0 mm microfiber filter. Then a subsample of 10 mL of this solution was passed through the FumoniTest column, which was washed with 10 mL of phosphate buffer solution followed by a second flow of 10 mL of phosphate buffer. The contents of the column were eluted with 1.0 mL of methanol (HPLC grade), collected, and mixed with 1 mL of developer. Finally, the fumonisin concentration in the grain was quantified using the FumoniTestTM and the Fluorometer VICAM, following the manufacturer’s protocols (Jorge et al., 2022).

Pedigree and genomic information were available for a total of 246 hybrids originated from 236 unique inbred lines (Table 1) belonging to two heterotic groups: 142 inbred lines from the Flint group, and 94 inbred lines from the Dent group (see Supplementary Figure S1 for more details about the heterotic groups). The inbred lines were divided into two groups of parents: P1, 30 inbred lines used as male; and P2, 211 inbred lines used as female. A few inbred lines (only 5) were used as both male or female; thus, these genotypes are found in both P1 and P2 groups.

For genotyping, the young leaves of 236 inbred lines were used for extracting genomic DNA via the CTAB method (Hexadecyltrimethylammonium bromide; Saghai-Maroof et al., 1984). The DNA quantification was done using a fluorometer, following the manufacturer’s instructions. The samples were shipped to the Genomic Diversity Facility of Cornell University (Ithaca, NY, United States) for genotyping-by-sequencing (GBS; Elshire et al., 2011). Using the Burrows-Wheeler alignment (BWA) tool (Li and Durbin, 2009), the sequences were aligned with the B73 reference genome (AGPv3). A total of 474,367 SNPs markers were available for analysis.

After applying quality control using TASSEL v.5.2.10 software (Bradbury et al., 2007), discarding those molecular markers with a Minor Allele Frequency (MAF) smaller than 10%, and a heterozygote’s proportion per locus above 10%, 73,083 polymorphic SNPs remained in the analyses. In addition, the SNPs with missing information were imputed using Beagle software (Browning and Browning, 2016). The SNP markers of the inbred lines were encoded as 0, 1, and 2, with 0 representing the allele with the minor frequency, one for the heterozygous, and two for the allele with the major frequency. The genotypes of the synthetic hybrids were built by combining the markers from the respective parental inbred lines using the expected value, which represents the mean allele dosage across parents for each marker.

In addition, there were also available phenotypic and pedigree information for an extra set of 127 hybrids totaling 373 when combined with the initial set of 246 hybrids that have both pedigree and genomic data (Table 1). This set of 127 hybrids was originated from 131 inbred lines belonging to three heterotic groups: 74 inbred lines from the Flint group, 49 inbred lines from the Dent group, and four inbred lines from the C group. The C group is intermediate to Flint and Dent, representing lines of several origins. It is well known there is a good combining ability between genotypes from the C group crossed with Dent or Flint testers (groups) (Silva et al., 2020). There was also information available regarding the pedigree of the ancestors of the inbred lines. Hence it was possible compute the pedigree-based relationship matrix for the hybrids ( A H ) and the corresponding matrices for the group of parents P1 ( A P 1 ), and P2 ( A P 2 ). It is worth noting that over 3 years, in the same experiment, there were hybrids whose parents were genotyped (246 hybrids) and hybrids whose parents were not genotyped (127 hybrids), due to limited resources of the breeding program.

The pedigree data was used to build the additive relationship matrix for P1 ( A P 1 ), P2 ( A P 2 ), and the hybrids ( A H ), according to Henderson (1976), which utilizes genealogical information to estimate genetic relatedness. This method is widely used in animal and plant breeding to quantify the genetic similarity among individuals based on their shared ancestry. The R package AGHmatrix (Amadeu et al., 2016) was used to build the additive relationship matrices for P1 ( A P 1 ), P2 ( A P 2 ), and the hybrids ( A H ). The genomic data was used to compute the realized genomic relationship matrix (or G matrix) for P1 ( G P 1 ), P2 ( G P 2 ), and the hybrids ( G H ) following VanRaden et al. (2008). VanRaden’s method estimates genetic relatedness based on actual genetic marker data, providing a more precise measure of genetic similarity by considering the genetic variation at specific loci across the genome. In addition, a) to having a better characterization of the matrix of genomic relationships, and b) to increasing the training set size by including non-genotyped inbred lines, c) the single-step relationship matrices were built for P1, P2, and the hybrid components to include information of the ancestry of the parents as in Khanna et al. (2022).

In general, the matrix that combines genomic and pedigree data is commonly named the H matrix (Velazco et al., 2019); however, since this research involves the prediction of hybrids, the use of H was reserved for denoting the genomic relationship matrix computed with the synthetic matrix of marker SNPs (mean of the allele dosage of the parents at each marker position). For this reason, the resulting matrix from the single-step procedure is named B-matrix. For example, the E + H model represents a linear predictor that includes the main effect of the environments and the main effect of the synthetic hybrid markers obtained as the mean across the marker information of the two parents involved in a specific crossing. This change was intended for an easier understanding for the reader. On the other hand, generically, the B (single step) matrix was built following Aguilar et al. (2010): B = w × A + 1 − w × G (1) where A is the pedigree relationship matrix, G is the genomic relationship matrix, and w is the weighting factor (i.e., the fraction of total additive variance not addressed by the markers; Velazco et al., 2019). Here, the B-matrix in Equation 1 can be computed considering different sources of information. For example, combining the H and the A matrices, or the G P 1 and the A P 1 matrices, the G P 2 and the A P 2 matrices, and their corresponding interactions with environments, etc. Further details about all the different combinations and model terms considered in this study are below.

A two-stage approach was considered to implement the prediction models. In the first stage, within-environment adjusted means (best linear unbiased estimation–BLUE) were obtained with the following model: y i k l = μ + h i + r k + b l k + ɛ i k l (2) where y i k l is the phenotypic value of the i ^th hybrid within the l ^th block of the k ^th replicate; μ is the constant, h i is the fixed effect of ith hybrid ( i = 1, 2, ., 373), r k is the fixed effect of replication k ( k = 1, 2), b l k is the random effect of the l ^th block within the k ^th replicate ( l = 1, 2, … ,16), with b l k   ∼   N 0 , σ b 2 , where σ b 2 is the block variance, and ɛ i k l is the residual random effect with ɛ i k l   ∼   N 0 , σ ɛ 2 , and σ ɛ 2 as the residual variance. A descriptive statistical analysis of the dataset was performed to assess the variance components and FUMO heritability for each environment and the results are presented in Supplementary Table S1. The analysis of the variance components employed a similar model to Equation 2 with the genetic effect of the hybrids considered as random such that h i   ∼   N 0 , σ h 2 , where σ h 2 represents the corresponding variance component. The first stage analysis was carried out using the statistical package ASReml-R v.4 (Butler et al., 2018) implemented in the statistical software R v.4.1.3 (R Core Team, 2022).

In the second stage, five different linear predictors were considered to model the trait response, and three different ways to account for relationships between hybrids using genomic (H, G _P1, G _P2), pedigree (A _P1, A _P2, A _P1P2), and the single-step (B, B _P1, B _P2, B _P1P2) information. Thirteen different prediction models resulted from combining the five linear predictors and the three different ways to introduce the relationships between hybrids (Figure 2; Supplementary Table S2). In principle, the first two linear predictors were implemented to model the hybrid performance directly, considering i) the main effects of the synthetic markers and ii) their corresponding interactions with environments. This model was also implemented to model the trait response using the pedigree matrix A of the hybrids and their corresponding interaction with environments. Hence, four different models can be constructed by combining the two linear predictors and the two ways to directly model the hybrids (M1, M2, M9, and M10 in Figure 2; Supplementary Table S2).

The other three linear predictors (3, 4, and 5) were considered to indirectly model the hybrids via the general and specific combining ability (GCA and SCA) terms. The GCA (P1+P2) term attempts to model the average effect of the parents involved in the crosses, while the SCA (P1 × P2) model term corresponds to the interaction of these effects. Similarly to the previous models, different sources of information were used to model the effects of the parents involved in the crosses: 1) molecular marker information ( G P 1 , G P 2 , and G P 1 P 2 ), 2) pedigree data ( A P 1 , A P 2 , and A P 1 P 2 ), and 3) the combination of these sources of information using the single-step method: B P 1 = w × A P 1 + 1 − w × G P 1 , B P 2 = w × A P 2 + 1 − w × G P 2 , and B P 1 P 2 = w × A P 1 P 2 + 1 − w × G P 1 P 2 . Further details are provided below to describe the election of the optimum value of the w hyper-parameter. In addition, the interactions of these terms with environments were also considered resulting in models M3, M4, M5, M6, M7, M9, M11, M12, and M13 (Figure 2; Supplementary Table S2).

Below, the linear predictors (1–5) and their corresponding components (Figure 2; Supplementary Table S2) are described.

For evaluating the model’s predictive ability, two cross-validation scenarios were implemented: 1) CV1, which consists of predicting hybrids that have not been tested at any of the observed environments, and 2) CV2, which consists of predicting the performance of hybrids that have been tested in some environments but not in others (incomplete field trials). Furthermore, two different ways of composing training sets represented in Figure 3 were considered. The traditional five-fold cross-validation scheme was implemented for the models that only use H, G, or the A matrices (models M1-M8; Figure 3A). This consists of a random five-fold (represented in Figure 3 by the hexagon), where four folds were used for model training, and the remaining fold was used as a testing set for evaluating the predictive ability. In this case, around 20% of the hybrids were predicted using 80% of the observed phenotypes. Ten replicates were considered to randomly assign phenotypes (CV2)/genotypes (CV1) to folds.

A slightly modified cross-validation approach was employed for the models incorporating the B-matrix (five models M9-M13), as depicted in Figure 3B. While maintaining the same five-fold partition, an optimization procedure was introduced to select the optimal w value to combining both covariance structures. In this setup, three folds were utilized for training the model, and the fourth fold was used to determine the optimum w value that yields the highest correlation between predicted and observed values. Regarding the fourth fold, a sequence of 21 values ranging from 0 to 1 (with increments of 0.05) for w was considered for constructing B = w × A + 1 − w × G matrices. After fitting the models using three folds, the correlation between predicted and observed values on the fourth fold was computed. The w value that returned the highest correlation was chosen to construct the B-matrix for predicting the fifth fold (testing set). Subsequently, the last fold was employed to assess the predictive performance of the models (Figure 3B). This approach allocated approximately 60% of the phenotypic data for model training, 20% for w optimization, and another 20% for evaluating the model’s predictive ability given the optimized w value obtained at the fourth fold.

The correlation between predicted and observed values was computed on a trial basis (within environments); thus, the choice of the w value can be determined using two approaches: 1) selecting the optimum value of w across environments and 2) selecting the optimum value of w for each environment. For the first case, the weighted average correlation across the three environments was computed following Tiezzi et al. (2017) (see Equation 8 for more details), then the w value returning the highest average correlation was used for predicting the fifth fold. For the second case, the w value that returned the highest correlation between predicted and observed values for each environment was selected for predicting the remaining fifth fold. Therefore, up to three different w values could result from selecting the optimum w value for each environment. This procedure was repeated 10 times for each one of the cross-validation scenarios (CV1 and CV2). The statistical analyses were performed in R statistical package version 4.2.1 (R core team 2022), using the package BGLR, version 1.1.0 (Pérez and De Los Campos, 2014).

The single-step model has shown improvements in predictive ability in plant breeding implementations (Legarra et al., 2009; Cappa et al., 2019; Oliveira et al., 2020). These improvements are attributed to the fact that i) combining genomic and pedigree information potentially results in a better matrix describing relationships between pairs of individuals, and ii) the ability to increase the training set size for model calibration (Oliveira et al., 2020). This study considered two scenarios for composing training sets to assess the factors contributing to improving predictive ability when considering the single-step model (B). The first scenario, named genomic only dataset (GOds), consists of only individuals with genomic and pedigree information (246 hybrids). While the second scenario, named pedigree and genomic dataset (PGds), increases the training set size by including phenotypic information of individuals with pedigree data (373 hybrids). Some of these have both genomic and pedigree information. For the two training set composition scenarios, in principle, the same training set partitions described in Figure 3 were considered; however, for the second scenario, the training set was augmented with the phenotypic information of non-genotyped inbreds for whose pedigree information was available to connect with genotyped inbreds.

The predictive ability was assessed on a trial basis by computing Pearson’s correlation between predicted and observed values within the same environment. The average correlation across environments was calculated by accounting for uncertainty and the sample size of the environments (Tiezzi et al., 2017) as follows in Equation 8: r ϕ = ∑ j = 1 J r j V r j ∑ j = 1 J 1 V r j (8) where r j is the Pearson’s correlation between the predicted and the observed values at the j ^th environment, V r j = 1 − r j 2 n j − 2 corresponds to the sampling variance, and n j is the corresponding number of observations. Furthermore, the across-environments mean squared prediction error (MSPE) was also computed as described in Equation 9. M S P E = ∑ j = 1 J ∑ i = 1 I y i j − y ^ i j T o t a l   p h e n o t y p e s (9) where y i j is the phenotypic value of the i ^th hybrid at the j ^th environment, and y ^ i j is the corresponding predicted value using the above described 13 models (M1-M13).

Two cross-validation scenarios (CV1 and CV2) were considered to assess the predictive ability of the prediction models. A total of 13 prediction models (M1-M13) were derived from the five linear predictors and the three different sources of information to modeling the different model terms. For models M1-M8, the training set was composed of four folds, while the fifth fold corresponds to the testing set. On the other hand, the models M9-M13 involve single-step parameterization. Thus, these require an optimization process to select the value of the w hyper-parameter that returns the highest predictive ability using three folds for model calibration and the fourth fold to evaluate different values for w (21 values ranging from zero to one in steps of 0.05). In both groups of models (M1-M8, and M9-M13), the same partitions were considered for a direct comparison.

In addition, for all models using only genotypes with genomic data (GOds) for composing the training sets, an alternative was considered by augmenting the four folds with phenotypic information of inbreds with pedigree data (PGds) only. Thus, in this case, it is possible to evaluate whether the improvements in predictive ability are because i) a better relationship matrix is computed when combining genomic and pedigree data, or ii) direct benefits of increasing the training set by including phenotypic information of non-genotyped inbreds, or iii) interaction of both events occurring at the same time.

Two different approaches were considered to find the optimum value of the w hyper-parameter used to blend/combine the A and the genomic (H or G) matrices according to the different linear predictors (1–5), the different ways for fitting the model terms (models M9-M13), and the different manners to compose calibration sets (i.e., GOds or PGds). The first case consists of optimizing w to return the highest weighted average correlation between predicted and observed values across environments, while the second approach focuses on finding the optimum w value for each environment such that up to three different values can be obtained. Figure 4 presents the weighted average (10 replicates) correlation across the three environments for both cross-validation schemes (CV1 and CV2), five prediction models M9-M13, two different ways to compose training sets (GOds and PGds), and both methods to find the optimum value of the w hyper-parameter (across and within environments).

As expected, when comparing these two approaches, for both prediction scenarios (CV1 and CV2) and both datasets (GOds and PGds), the predictive ability always improved when selecting w within environments for models M9-M13. Regarding the training set compositions, considering only information of genotyped individuals (GOds), under the CV1 cross-validation scheme, the average predictive ability ranged from 0.11 to 0.21, while between 0.09 and 0.18 for CV2 when performing the selection of the w value across environments (Figure 4A). On the other hand, when the w value was optimized within environments, the corresponding average predictive abilities increased, ranging from 0.18 to 0.26 (CV1) and 0.14 to 0.21 in (CV2; Figure 4A). The model that showed the highest predictive ability, for the two approaches (within and across environments), considering the optimization of the w hyper-parameter in both prediction scenarios (CV1 and CV2) and considering training sets of only genotyped individuals (GOds) was M10 (Figure 4A).

Augmenting the training set considering phenotypic information of individuals with pedigree data (PGds) when using the within environments w optimum value, the average predictive ability ranged from 0.20 to 0.32 and 0.12 to 0.21 for CV1 and CV2, respectively (Figure 4B). The model that showed the highest predictive ability for CV1 using w optimum across environments was M13. Nonetheless, when using the within environments w optimal, the highest predictive ability was shown with M10 (Figure 4B). These results indicate that the manner of computing the optimized B-matrix influences the election of the best prediction model. For CV2, the model that showed the best predictive ability, for both approaches (GOds, and PGds), within and across environments optimization, was M13. Due to the results presented, the w optimum value computed within environments was used on the models considering the B matrices (M9-M13) for comparing the predictive ability of the 13 models (M1-M13) for GOds and PGds training sets (see Figure 4).

For an easy visualization, Figure 5 presents the heatmap of the within and across-environments predictive ability and the across-environments mean squared prediction error (MSPE) of the 13 models (M1-M13) for the two predicting scenarios (CV1 and CV2) and the two manners to compose training sets (GOds and PGds). As mentioned before, for models M9-M13, only the results derived from the within environments optimization are considered.

The predictive abilities varied considerably through environments for both prediction scenarios (CV1 and CV2) and both datasets used for model training (GOds, and PGds). Overall, there was an improvement in the across-environments weighted predictive ability and a reduction of the MSPE when the B-matrix was considered (higher correlations and reduced MSPE) in models M9-M13. Note that despite using phenotypic information from genotyped individuals (GOds) or combining it with also individuals with only pedigree data for increasing training set size (PGds), the models including the B-matrix improved the predictive ability compared with those that only used the genomic information (H or G matrices) in M1-M5 models. In general, the models presenting the best results, within and across environments, are those that considered the interaction between the genetic (H, G, A, or B) and the environmental component (E).

In the previous section, across environments, mixed results were shown with the increasing of the training set by adding those individuals with pedigree information only. For CV1, the best results were observed when the training set was augmented (0.321) via the PGds data set in comparison with the GOds set (0.280). Recalling that in both cases, for the most promising models (M10-M13), the genomic information and the pedigree data are combined via the B matrices (after an optimization process). However, the only difference is that for the GOds the training set is composed for only those individuals that have both genomic and pedigree data, while for the PGds the previous data set is augmented by adding the phenotypic information from individuals with only pedigree information.

For CV2, an opposite trend was observed, with GOds returning better results (0.229) than PGds (0.171) when considering the most promising model (M10) from the previous cross-validation scenario. However, in this case, for PGds M13 returned the highest correlation (0.233). In any case, for M11-M13 models, slightly improvements were observed when using PGds vs. GOds. These differences might not be significant; thus, no clear advantages were observed with the increasing training set size as for the CV1 case. Figure 6 depics a direct comparison of the effects in predictive ability i) from combining genomic and pedigree information (GOds), ii) with also an increased training set size (PGds) with respect to the models based on only a single data type (H, G or A used to describe similarities between pairs of individuals).

More specifically, in Figure 6, the results of the within and across environments predictive ability of the five linear predictors are depicted for single H (M1-M2), G (M3-M5) and combined with A through B (M9-M13) covariance structures considering phenotypic information of only individuals with genomic and pedigree data (GOds), and adding information from individuals with only pedigree data (PGds).

Under the CV1 scheme, for the linear predictor 1, on average moderate improvements are observed by combining the matrices H and A (especially for E2) compared with only using the H matrix to model relationships between pairs of individuals (M9 vs. M1). Also, in average, slight improvements were observed when the training set size was increased with individuals with only pedigree data (PGds) compared to the case where only the phenotypic information of individuals with both genomic and pedigree data (GOds) was used for model calibration.

A similar pattern was observed for the second linear predictor, where the interaction with environmental factors was included using the same covariance structures (H vs. B) contrasting models M2 and M10. In addition, the linear predictors 3-5 showed similar patterns. In all the cases, very low correlation values were observed for E1, suggesting this environment’s high influence contributing to the G × E interaction.

Regarding the CV2, more pronounced improvements in predictive ability were found by combining H and A matrices for linear predictors 1 and 2. However, these decreased when the information from non-genotyped was also included in the calibration sets, except for E2, where it practically remained the same. For the linear predictors 3 and 4, some improvements were observed with the PGds compared to GOds, while with the most complex linear predictor (5), mixed results were found depending on the environment.

This study introduces an approach for optimizing the hyper-parameter w using a novel method. To find the optimum w value, a cross-validation strategy using four of the original five folds was implemented (three folds for model fitting and one-fold to optimize w for 21 different values evenly spaced between 0 and 1 in steps of 0.05). Then the phenotypic information of the fifth fold was predicted considering the obtained w optimum value. This strategy offers a robust method for selecting the hyper-parameter w, returning satisfactory results regarding the accuracy of genomic predictions, which is particularly crucial for complex traits such as FUMO. This suggests that using the fourth fold for B-matrix optimization can be a good strategy to guide the election of the w hyper-parameter.

Available studies have shown that the election of the w optimum value varies according to the different characteristics being evaluated (Liu et al., 2011; Ashraf et al., 2016; Velazco et al., 2019). Nonetheless, our study presents a significant contribution by demonstrating the importance of selecting the optimal w by environment in unbalanced datasets, an approach not yet utilized in previous research. Until now, the studies have focused on selecting the w optimal value across environments (i.e., a common optimum value for all the environments). However, individuals have different genetic responses when changing environmental conditions, with non-static variation (e.g., temperature, precipitation, humidity) being the main source of genotype-environment interaction for maize (Cullis et al., 2000; Kleinknecht et al., 2013; Dias et al., 2020; Krause et al., 2022).

Our results reveal that, particularly in situations with unbalanced datasets, the difference in gene expressions of individuals between environments can significantly influence the election of the w optimal weighting factor. This finding has important implications for multi-environment and multi-trait studies, suggesting that researchers should consider not only the selection of w by trait, as is customary, but also by environment. The change in the value of w with unbalanced data occurs due to the alteration in the relative contribution of genomic and pedigree information to the genetic structure, which may necessitate a reevaluation of the hyper-parameter w to ensure that the modeling of genetic relationships is accurate and appropriate for the specific data in question. Consequently, our study demonstrates that the best strategy for building the B-matrix should be optimizing the w hyper-parameter within in environments, which is particularly relevant for crop improvement in tropical environments where G × E interactions are pronounced, and datasets are often unbalanced.

Over the decades, breeders have sought strategies to reduce the FUMO content in grain (Duvick, 2001; Eller et al., 2008; Pádua et al., 2016; Santiago et al., 2020; Butoto et al., 2022). The USFDA Center for Food Safety and Applied Nutrition (USFDA, 2021) recommends a maximum of 2–4 mg kg^-1 of contamination of FUMO in corn products. In Brazil, the contamination limit is 2 mg kg^-1 of these toxins in grains (Agência Nacional de Vigilância, 2011). Exceeding these limits prevents the exportation and national marketing of maize batches, leading to substantial economic losses for producers. Therefore, in addition to high grain yield performance, developing cultivars resistant to FUMO is a primary objective. However, phenotyping this trait is laborious and expensive (Bush et al., 2004). Consequently, implementing genomic selection becomes an important alternative strategy for developing fumonisin-resistant genotypes. Comparing genomic to traditional phenotypic selection for FUMO Butoto et al. (2022) found that both methods performed similarly. Nevertheless, the authors highlight that genomic selection has the potential to be more efficient than phenotypic selection due to the employment of cheaper and faster genotyping methods.

To our knowledge, no studies in the literature have combined genomic and pedigree information to construct relationship covariance structures for predicting FUMO. Our results demonstrate that the use of the B-matrix improved predictive ability in all tested linear predictors for both prediction scenarios and datasets (GPds and POds) compared to conventional GS models. This improvement in predictive ability is crucial for complex traits like fumonisin resistance, where accurate predictions can significantly impact breeding program efficiency. As observed, the use of a single-step approach was superior to the use of only the H matrix, regardless of whether the training set size was increased with un-genotyped individuals. Since some markers may not be in linkage disequilibrium with QTLs, when combining the A matrix and H or G matrix in a relationship the pedigree information may have contributed to capturing associations between causative alleles due to common ancestral identity, improving predictions models (Velazco et al., 2019). This finding has important practical implications, especially in situations where genotyping costs are a constraint. The B-matrix approach allows increase the size of the training set by including un-genotyped individuals through pedigree information, potentially leading to more robust and accurate predictions.

The observed improvements in prediction accuracy using the B-matrix align with findings from previous studies by Imai et al. (2019) and Cappa et al. (2019). These authors found that in terms of prediction accuracy, the B-matrix matched or surpassed the use of only the G matrix or the A matrix only. Additionally, over the years, several other authors have reported enhanced prediction accuracy with the use of the B-matrix in plant and animal species (Ashraf et al., 2016; Gao et al., 2012; Koivula et al., 2012; Lourenco et al., 2015; Pérez-Rodríguez et al., 2017). Furthermore, the single-step method allows the utilization of larger phenotypic datasets compared to the GBLUP method, as demonstrated herein for maize. Moreover, all the tested models corresponding to the different linear predictors reached the lowest MSPEs when using the B-matrix in comparison to those using only one data type to model similarities between pairs of individuals (H, G, or A). These results are consistent with those obtained by Velazco et al. (2019). The selection of prediction models based on minimizing MSPE has been recommended because this statistic considers both the precision and bias of the models (Vitezica et al., 2011; González-Recio et al., 2014).

Mainly in animal breeding, the use of the B-matrix in genomic prediction has been widely discussed (Martini et al., 2018; Macedo et al., 2020; Mäntysaari et al., 2020; Masuda et al., 2021). However, this methodology has not yet become very popular in plant breeding (Oliveira et al., 2020). Our results demonstrate the advantage of using the single-step approach in the intermediate stage of a breeding program in two different contexts. The first context is when the program does not have all individuals genotyped (PGds), for example, due to new materials inserted in the pipeline. In this case, the B-matrix enables the construction of the relationship matrix incorporating all individuals improving the predictive capacity of the models. The simultaneous use of genotyped and un-genotyped individuals relies on projecting genomic relationships to un-genotyped individuals based on the conditional distribution of breeding values for un-genotyped and genotyped individuals (Legarra et al., 2009). The genomic relationships can improve the pedigree relationships, while the un-genotyped individuals provide more phenotypic information. As a strategy, de Oliveira et al. (2020) used the B-matrix in multi-trait multi-environment genomic prediction models due to the lack of genotypic information for some evaluated maize hybrids. Likewise, Cappa et al. (2019) noted that the utilization of the single-step method resulted in elevated prediction accuracies and reduced bias of the genetic component of unobserved (non-phenotyped) but genotyped individuals compared to the standard GBLUP by using additional phenotypic information from non-genotyped individuals.

In our study, including phenotypes from un-genotyped individuals increased the predictive ability in most but not all the tested scenarios. In agreement with our findings, in a study with a broiler population, Hidalgo et al. (2021) found that the two most recent years of pedigree, phenotypic, and genomic data were sufficient to maintain prediction accuracies in selection candidates (i.e., the last generation of individuals), adding phenotypes of un-genotyped individuals from previous years did not increase the accuracy. Lourenco et al. (2014) stated that distant ancestors have minor contributions, explaining the null or marginal increase in predictive ability, and sometimes their inclusion can deteriorate predictive ability.

The second context consists of companies with genomic and pedigree information for all individuals (GOds). Our results showed that by combining both sources of information for modeling covariance structures, the predictive ability of the models increases, helping to prevent (discard) of advancing to the next stages of the program the most susceptible hybrids to FUMO. The increased predictive ability can be explained because the markers do not capture all the genetic variance; blending the genomic relationship matrix with a portion of the pedigree relationship matrix (with an optimal value w ) implicitly fits a residual polygenic effect in the statistical model via the modified genomic relationship matrix. Although the genomic relationship matrix accounts for most of the genetic variation, its combination with the pedigree relationship matrix increases the captured genetic variance. Furthermore, the combined use of pedigree and genomic information in the B-matrix can help mitigate the limitations of using the G matrix alone, such as incomplete linkage disequilibrium between markers and causal variants, especially in cases where the trait is influenced by rare or less common genetic variants, or imperfect genomic data. By leveraging the complementary information from both sources, the B-matrix provides a more comprehensive and accurate representation of the genetic relatedness among individuals, ultimately enhancing the predictive ability of the models. Other researchers have also observed that combining genomic and pedigree information optimizes genomic prediction for complex traits (Crossa et al., 2013; Basnet et al., 2019; Velazco et al., 2019).

Our study employed two cross-validation schemes, CV1 and CV2, to evaluate the predictive ability of various models. The CV1 scheme, which predicts untested hybrids in evaluated environments, generally showed higher predictive abilities compared to the CV2 scheme. This finding was contrary to our initial expectations based on similar studies in the literature (Jarquin et al., 2021; Khanna et al., 2022; Persa et al., 2021). One of the factors that could contribute to this unexpected outcome is the unbalanced nature of our dataset that could potentially impact the effectiveness of the CV2 approach, which relies on information from other environments. These findings underscore the need for further research to understand the factors influencing the relative performance of CV1 and CV2 schemes in different contexts, particularly for traits showing strong environmental influences. Finally, the breeding programs are also interested in developing superior cultivars that respond favorably to diverse environmental conditions (Jarquin et al., 2021). Our results were consistent with those obtained by several researchers that observed better predictive ability in models that consider genetic and environmental interaction effects (Jarquin et al., 2014; 2021; Lado et al., 2016; Basnet et al., 2019; Khanna et al., 2022), further emphasizing the importance of accounting for G × E interactions for reaching improved prediction accuracies of complex traits in tropical maize breeding programs.