Genomic Selection Using Pedigree and Marker-by-Environment Interaction for Barley Seed Quality Traits From Two Commercial Breeding Programs

With the current advances in the development of low-cost high-density array-based DNA marker technologies, cereal breeding programs are increasingly relying on genomic selection as a tool to accelerate the rate of genetic gain in seed quality traits. Different sources of genetic information are being explored, with the most prevalent being combined additive information from marker and pedigree-based data, and their interaction with the environment. In this, there has been mixed evidence on the performance of use of these data. This study undertook an extensive analysis of 907 elite winter barley (Hordeum vulgare L.) lines across multiple environments from two breeding companies. Six genomic prediction models were evaluated to demonstrate the effect of using pedigree and marker information individually and in combination, as well their interactions with the environment. Each model was evaluated using three cross-validation schemes that allows the prediction of newly developed lines (lines that have not been evaluated in any environment), prediction of new or unobserved years, and prediction of newly developed lines in unobserved years. The results showed that the best prediction model depends on the cross-validation scheme employed. In predicting newly developed lines in known environments, marker information had no advantage to pedigree information. Predictions in this scenario also benefited from including genotype-by-environment interaction in the models. However, when predicting lines and years not observed previously, marker information was superior to pedigree data. Nonetheless, such scenarios did not benefit from the addition of genotype-by-environment interaction. A combination of pedigree-based and marker-based information produced a similar or only marginal improvement in prediction ability. It was also discovered that combining populations from the different breeding programs to increase training population size had no advantage in prediction.

With the current advances in the development of low-cost high-density array-based DNA marker technologies, cereal breeding programs are increasingly relying on genomic selection as a tool to accelerate the rate of genetic gain in seed quality traits. Different sources of genetic information are being explored, with the most prevalent being combined additive information from marker and pedigree-based data, and their interaction with the environment. In this, there has been mixed evidence on the performance of use of these data. This study undertook an extensive analysis of 907 elite winter barley (Hordeum vulgare L.) lines across multiple environments from two breeding companies. Six genomic prediction models were evaluated to demonstrate the effect of using pedigree and marker information individually and in combination, as well their interactions with the environment. Each model was evaluated using three cross-validation schemes that allows the prediction of newly developed lines (lines that have not been evaluated in any environment), prediction of new or unobserved years, and prediction of newly developed lines in unobserved years. The results showed that the best prediction model depends on the cross-validation scheme employed. In predicting newly developed lines in known environments, marker information had no advantage to pedigree information. Predictions in this scenario also benefited from including genotype-by-environment interaction in the models. However, when predicting lines and years not observed previously, marker information was superior to pedigree data. Nonetheless, such scenarios did not benefit from the addition of genotype-by-environment interaction. A combination of pedigree-based and markerbased information produced a similar or only marginal improvement in prediction ability. It was also discovered that combining populations from the different breeding programs to increase training population size had no advantage in prediction.

INTRODUCTION
In the last two decades, plant breeding programs have advanced from conventional approaches involving visual selection and trait screening over several generations to potentially faster methods by means of marker technologies . Genomic selection, the concept of using dense molecular genome-wide data to predict the performance of individuals in a breeding population (Meuwissen et al., 2001), has become more popular in cereal breeding in recent years. Empirical studies in barley (Zhong et al., 2009;Sallam et al., 2015;Nielsen et al., 2016;Thorwarth et al., 2017), wheat (Triticum aestivum; Crossa et al., 2010;Dong et al., 2018;Haile et al., 2018;Norman et al., 2018), maize (Zea mays; Zhao et al., 2012;Shikha et al., 2017;Vélez-Torres et al., 2018), rice (Oryza sativa; Spindel et al., 2015;Xu et al., 2018), and sorghum (Sorghum bicolor;De Oliveira et al., 2018;Hunt et al., 2018) have recently all proven that with current advances in the development of high-density array-based DNA marker technologies and reduced costs, genomic selection has become an important tool in cereal breeding.
Prediction models that use realized relationships based on marker information lead to a substantial increase in the prediction accuracies of complex traits compared to those using relationships based on pedigree information (VanRaden et al., 2009), and this has been observed in several genomic selection studies (Crossa et al., 2010;Albrecht et al., 2011;Burgueño et al., 2012). It is assumed that by using the markerbased related matrix, it is possible to account for Mendelian segregation of alleles along with full-sibs, which are more or less expected due to random chance. Hence, while pedigree-based breeding values of unphenotyped lines will only reflect midparent genetic contributions, genomic-estimated breeding values of unphenotyped full sibs will reflect genetic differences caused by Mendelian sampling (Velazco et al., 2019).
Other researchers have found less conclusive results, which show that in some instances, pedigree-based relationships may perform similarly to or even better than marker-based relationships in terms of prediction accuracy (Juliana et al., 2017;Hunt et al., 2018). This is particularly true when the pedigree data are very precise and include several generations (Juliana et al., 2017). Numerous studies have also reported the benefits of jointly using marker and pedigree-based relationships in prediction models (Crossa et al., 2010;Albrecht et al., 2011;Burgueño et al., 2012;Rodríguez-Ramilo et al., 2014;Sukumaran et al., 2017). The inclusion of both marker and pedigree has the potential to increase prediction performance because the pedigree information can account for variations that may not be captured by markers at population and family level. Hence they provide better estimates of genetic variations, thereby improving predictive performance (Velazco et al., 2019). However, it has also been demonstrated that when high-quality genome-wide marker and pedigree data are available, it is not always expected that the combination of marker and pedigree information will outperform the use of marker information alone in mixed models due to information matrix redundancy (Albrecht et al., 2011).
In the past, genomic prediction models have been developed to suit single environment predictions. However, in most plant breeding programs, multi-environment trials play a crucial role in the assessment of lines to determine the performance of stability across environments of the prospect variety. The effects of genes on traits are usually influenced by environmental conditions, leading to high levels of genotype-by-environment interactions (Burgueño et al., 2012). Genomic selection models have therefore been extended to fit multi-environment settings. The incorporation of genotype-by-environment interactions (marker-by-environment; GxE and pedigreeby-environment; AxE) in prediction models has been found to increase the performance of such models in comparison with single environment models (Burgueño et al., 2012). Furthermore, it has been shown that prediction accuracy in barley increases when considering marker and marker-byenvironment interactions using ridge-regression best linear unbiased prediction to extend genomic selection to multiple environments (Oakey et al., 2016). However, when high dimensional genetic and environmental variables are involved, modeling all possible interactions becomes a challenge. Other studies (Jarquín et al., 2014;Pérez-Rodríguez et al., 2015) have also shown that incorporating interactions between genetic information and environment increased prediction accuracy by a considerable margin. In similar studies, interactions between genetic information based on both pedigree and marker and on 18 environment variables increased the prediction ability (Sukumaran et al., 2017).
The present study assessed the performance of prediction models on four seed quality traits in winter barley from two breeding companies in Denmark: Nordic Seed (hereafter referred to as NS) and Sejet Plant Breeding (hereafter referred to as SJ). Specifically, this study: (i) evaluated the individual and combined effect of marker and pedigree data on prediction ability under different cross-validation schemes, (ii) determined genetic (marker and pedigree) and environment-interaction effects on prediction ability under different cross validation schemes, and (iii) assessed the advantage of increasing training set size by combining lines from different breeding programs.

Phenotypic Data
Data were provided by two commercial breeding companies in Denmark: NS and SJ (Supplementary Material Table 1: "Phenotypic data"). Line records (1,850 in all) were obtained from 484 commercial elite two-row winter barley lines from NS and evaluated over 3 years (2015, 2016, and 2017) across three locations in Denmark. All three locations were the same for the 3 years. Data from SJ consisted of 964 records from 428 lines and phenotyped for protein content and test weight. Seed fraction weight records for SJ were measured in the Department of Plant and Environmental Sciences at University of Copenhagen, Denmark. The 428 lines from SJ were also grown in 3 years (2015, 2016, and 2017) across two locations in Denmark (except in 2016, where lines were evaluated at only one location).
At each location in the breeding companies, lines were grown in individual experiments consisting of 5 to 20 small trials with the number of lines per each trial ranging between 25 to 30. Two to 6 lines were used as checks in each of the experiment. Each trial was a randomized block-design with three replicates per each line. Management of experiments were the same for all locations within the specific breeding companies. Due to this, there is the possibility of differences in the management practices observed in the different breeding programs.
Five lines were common between the two companies. In all, 2,814 records from 907 lines were used in the analysis. On average, there were 14 distinct environments (companyyear-locations). The number of observations recorded in each environment is presented in Table 1. Each line record was phenotyped for protein content (% dry matter), test weight (kg/hL), percentage weight of seed fractions above 2.5 mm in diameter (hereafter referred to as SF_abv2.5), and seed fractions above 2.2 mm in diameter (hereafter referred to as SF_abv2.2). Protein content and test weight were measured using a FOSS Grain Analyzer instrument based on near-infrared transmission technology. The seed fractions were obtained using a SORTIMAT screening instrument to separate seed samples of approximately 100 g into four seed size classes of >2.8, 2.8-2.5, 2.5-2.2, and <2.2 mm in diameter. The seed size classes were divided into two main classes: those above 2.5 mm in diameter (>2.8 and 2.8-2.5 mm classes) and those above 2.2 mm in diameter (>2.8, 2.8-2.5, and 2.5-2.2 classes). Their percentage weights were calculated and used in the analysis.

Genotype information
Lines from both companies were genotyped with the 9K iSelect single-nucleotide polymorphism (SNP) chip by TraitGenetics. Filtering out monomorphic sites and more than 10% missing markers reduced the SNP markers to 4,830 for analysis (Supplementary Material Table 1: "Genotypic data"). The remaining missing SNPs were imputed with "synbreed" using the "random" algorithm, where missing values for a marker were sampled from the marginal allele distribution of that marker (Wimmer et al., 2012). The realized genomic relationship matrix based on marker information was computed (VanRaden, 2008).

Pedigree Information
Pedigree information on all the lines was provided for the analysis. Pedigree information of the 484 lines in the NS data extended to five generations. It comprised 80 subfamilies in total, with progenies per subfamily (full siblings) ranging between 2 and 24, and 19 lines representing individual crosses (Supplementary Material Table 1: "Pedigree matrix_NS data"). The pedigree information for all 428 lines from SJ extending up to seven generations comprised 106 subfamilies in total (Supplementary Material Table 1: "Pedigree matrix_SJ data"). The number of progenies per subfamily ranged from 2 to 13 full siblings, with 67 lines with individual crosses (no siblings). This shows that although the number of subfamilies within NS was small compared to SJ, the subfamily units within NS had more siblings than SJ.

Statistical Models
Data were analyzed using mixed linear models, where the main effect of lines, environments, markers, pedigree, and their interactions with environments were modeled using random covariance structures. Since data on line records, marker and pedigree information were used, the line was included as main effect in all models so as to partition the variance components into additive (marker and pedigree), and non-additive (residual genetic) effects. The seed quality traits analyzed were protein content, test weight and two seed fraction weights. The models used for evaluation are briefly described below.

Model 1: ELG
In this model, a linear mixed model that accounts for main environmental effects, line effects, and marker effects was fitted as expressed below: where y ij , is the phenotype, µ is the intercept, E i is the random effect of the ith environment (company-year-location), L j is the random effect of the jth line, and ε ij is the residual error. g j is the random "genomic breeding value" of the jth line. This accounts for the marker-based relationships of the lines by following a multivariate normal density, with zero mean and covariance matrix Cov g = Gσ 2 g , where G is the realized genomic relationship matrix derived by the VanRaden formula XX 2 p m (1−p m ) (VanRaden, 2008), where X is the centered and standardized marker matrix, p m is the allele frequency of the mth marker, and σ 2 g is the marker genetic variance. All other random effects in the above model follow normal independent and identical distributions, such that E i ∼ N(0, Iσ 2 E ), L j ∼ N(0, Iσ 2 L ), and ε ij ∼ N(0, Iσ 2 ε ).

Model 2: ELG-GxE
This model was obtained by extending model 1 by incorporating the interaction between the random effects of the markers and the environments, such that the model accounts for the main effects of the environments, the main effects of the lines, the main effects of the markers, and the interactions between the markers and the environments. The interaction term is incorporated as a Hadamard product of two covariance structures (E i and g j ) describing relationships between the environments and the lines based on the marker information. Model 2 is expressed as: where gE ij is the interaction between the jth line based on the marker data in the ith environment. The other components of the model are as defined previously. The random vector containing the interaction term is assumed to follow a multivariate normal distribution gE ij ∼ N(0, Z g GZ g • Z E Z E σ 2 gE ), where Z g is the incidence matrix for the vector of marker effects, Z E is the incidence matrix of the environmental effects, σ 2 gE is the variance of the interaction term gE ij , and • represents the cell-by-cell or Hadamard product between the two matrices.
Model 3: ELA Model 3 was obtained by replacing the random effect of the marker covariate in model 1 with the random additive relationship of the lines based on the pedigree information (a). Model 3 is represented as: where a j is the random additive effect of the jth line, with a covariance structure based on the pedigree relationships between the lines, so that a j follows a multivariate normal density with zero mean and covariance matrix Cov (a) = Aσ 2 a , where A is the additive pedigree-relationship matrix, and σ 2 a is the pedigree genetic variance.

Model 4: ELA-AxE
Model 4 was obtained by adding the interaction between the lines based on their pedigree information (a) and the environments (E) to model 3. As in model 2, the interaction was incorporated as the Hadamard product of the covariance structures for E i and a j and is thereby represented as: where aE ij is the interaction term between the jth line based on the pedigree information in the ith environment assumed to follow a multivariate normal distribution of the form . Z a is the incidence matrix for the vector of pedigree effects, and σ 2 aE is the variance of the interaction term aE ij .
Model 5: ELGA This model was obtained by combining the components of model 1 and model 3, such that the model accounts for random effects of the environments (E), lines (L), marker (g), and pedigree (a). Model 5 is represented as: where all the terms of this model are as previously described.

Model 6: ELGA-GxE-AxE
Model 6 is literally the combination of the components in models 1 to 5, making it the most comprehensive model. With all the terms previously defined, the model is expressed as.
All the above models were fitted in R using the Bayesian Generalized Linear Regression (BGLR) package (Pérez and de los Campos, 2014; de los Campos and Pérez-Rodríguez, 2018).

Assessing Prediction Performance for Different Cross-Validation Strategies
Genomic prediction models were initially fitted with data from each company and a combined dataset to estimate the variance components for each model. To assess the performance of all the models for each breeding program separately, three crossvalidation schemes were implemented.
(i) CV1: Random 10-fold cross validation. This scheme mimics breeding problem of predicting the performance of newly developed lines that have not yet been observed in any environment . Here, lines were randomly partitioned into 10 folds and nine folds (90% of the lines) were used to train the model to predict the remaining one fold (10% of the lines). The partitioning was done such that all line records from the same genotype appear in the same fold. The process was repeated ten times until each fold had been used once as the validation set. The correlation between the observed phenotypic measures and the estimated breeding values (k) for each genotype (r y,k ) was estimated and the average of the correlation coefficients from the ten folds reported as the prediction ability. The estimated breeding valuek is defined in respect to the specific model used. Standard error from the mean of the 10-folds correlations were recorded. (ii) CV0: Predicting unobserved or new year. This scheme allows the prediction of crop performance of all lines in unobserved (future) year. Also referred to as forward prediction or leave-one-year-out, this scheme allows the prediction of environments in a given year using all the environments from the previous years. In this study, data from 2015 and 2016 were combined to form the training set, and year 2017 was used as the validation set. For prediction with SJ data, there were 303 and 125 in the training and validation sets, respectively, while for NS, there were 317 and 167 lines in the training and validation sets, respectively. This scheme provides a more realistic problem often faced by breeders which has to do with unobserved years. The success of this strategy depends not only on the number of related lines in the training and validation sets but mostly on the correlation between the environmental conditions of the unobserved years and the previous years . (iii) CV00: Predicting newly developed lines in new year. This scheme aims at predicting the performance of all new lines in the future year. This scheme is similar to CV0, however, the lines to be predicted have never been observed in the years previously. This scheme applies to breeding situations where most of the materials have not been observed in previous fields but their performance needs to be estimated for the next year Howard et al., 2019). Here lines that were evaluated only in 2017 but not in any of the previous years were used as validation set while the rest were used as the training set. For SJ predictions, there were 348 and 80 lines in the training and validation sets, respectively, while for NS there were 374 and 110 lines in the training and validation sets, respectively.
In CV0 and CV00, the procedure did not involve random partitioning, and hence implemented once. The prediction ability was estimated as correlation between the predicted values and the observed phenotypic values. Predictions using individual breeding program data were referred to as withincompany predictions.
Subsequently, we assessed the effect of combining different populations originating from different programs to increase training population size in prediction. We referred to these scenarios as across-company predictions while those using data specific to the different breeding programs are referred to as within-company predictions. In the across-company predictions, data from both breeding programs were combined and lines were randomly divided into ten folds. Folds were created such that nine of them contain lines from both companies, while the remaining fold contain either (i) only lines from SJ as the validation set or (ii) only lines from NS as the validation set. All models previously discussed were modified by including breeding program as a fixed effect in the across-company predictions. Predictions derived from this approach were compared to the ones derived in the CV1 of within-company predictions discussed above.

Estimating Heritability From Unbalanced Data
Broad-sense (H 2 ) and narrow-sense (h 2 ) heritabilities were estimated, respectively, as: where phenotypic variance was computed as, σ 2 For models containing marker relationships alone, the value σ 2 a = 0 was used, and for models containing pedigree relationships alone, the value σ 2 g = 0 was used. Furthermore, the values σ 2 gE = 0 and σ 2 aE = 0 were used for models without environmental interactions. Given that the analysis was based on unbalanced data, the equation was followed with the definition m h = n n i=1 where n is the number of genotypes, m i is the number of environments for each genotype, and p i is the number of plots across environments for each genotype (Holland et al., 2003).

Descriptive Statistics
The lines used in this study came from active barley breeding programs run by two breeding companies. In each breeding company, lines were grown between one and three locations over three years, comprising a total of 14 environments (companyyear-location combinations), with NS having nine environments and SJ five environments. Table 2 gives the summary description of the respective traits. Overall, the average protein content was 10.48% (SD = 0.83), and that of NS and SJ was 10.54% (SD = 0.69) and 10.36% (SD = 1.04), respectively. The average and standard deviation estimates of the remaining traits given in the table shows that seed fraction SF_abv2.5 was the most variable, followed by protein, and with the least variable being seed fraction SF_abv2.2. Figure 1 presents a boxplot of all the traits analyzed for each of the 14 environments in the study. Except for protein content, all the traits showed higher values in 2015 than in 2017 for both breeding companies. For protein and test weight, the variations in the environments seem dependent on both year and location, while for the seed fractions it seems more dependent on year in both breeding companies. The boxplot shows significant environmental effects in all the traits. Analysis of variance of the 14 environments shows F values of 510, 443, 260, and 163 for protein, test weight, SF_abv2.5, and SF_abv2.2, respectively.

Relationship Between Lines
The degree of relatedness of the lines based on the pedigreebased and marker-based information are depicted in the heat maps in Figure 2. The values of the pedigree-based relationship matrix are derived from the probability of identity by descent (IBD) and therefore the values are between zero and one. Values of zero mean no relationship and one indicate an identical relationship. Unlike the pedigree-based matrix, values of the marker matrix are composed of negative and positive values. The negative relationships are explained from the centering of the marker covariates, leading to centering of the entire marker-based matrix such that the sum of all elements in the matrix is zero. Negative values in the marker-based relationship matrix mean that detection of an allele in one line makes it less likely to be detected in the other line, zero means a lack of dependence and positive values indicate an increased likelihood of the allele being detected in the other line.
The pedigree heat map shows fewer groupings due to the assumptions of IBD that allow full siblings derived from the same cross to be similar. The marker heat map, however, shows many more groupings. This was to be expected because the marker matrix accounts for Mendelian segregation of alleles, such that full siblings that are more or less expected due to random sampling can be distinguished. The heat maps of both matrices showed a closer relationship between NS lines (Figure 2A) than SJ lines (Figure 2B), as is evident from the greater density of blue and yellow shades in the NS heat maps compared to SJ. Moreover, the mass of SJ lines at zero (lack of dependence), as shown in the color key density plot for the marker matrix, were higher than that of NS, indicating a higher relatedness of NS lines compared to SJ lines.

Variance Component Estimates
A summary of the variance components for the studied traits estimated from all six models using combined data and individual company data are presented in Tables 3A-C, respectively. The main effects of environments consistently explained most of the total variance for all traits. Using SJ data alone, environmental variance (E) explained between 42% (in SF_abv2.2 ELGA-GxE-AxE model) and 81% (in protein ELG model) approximately of the total variance (Table 3B). Environmental variance explained between 36% (in SF_abv2.2 for ELA-AxE model) and 63% (in protein ELG model) approximately of the total variance using NS data (Table 3C). Due to the large proportions of variance explained by the environments in the different models, the proportional contribution of each of the random effects relative to the total phenotypic (within-environment) variance was estimated, as suggested in literature (Jarquín et al., 2014;Pérez-Rodríguez et al., 2015).
Depending on the model and trait, lines explained between 8% and 19% (with SJ data), and between 8.5 and 28% (with NS data) of the within-environment variance. Marker information explained between 12.9 and 34.7% of the within-environment variance, while pedigree information explained between 12.8 and 44.5% (for SJ data), and depending on the trait and the inclusion (or not) of interaction effects. With NS data on the other hand, marker information explained between 15.4 and 43.8% while pedigree information explained between 13 and 49.5% of the within-environment variance.
There was a decrease in the estimated within-environment error variance (R) when interaction terms were included in the main effect models for all traits. Considering protein in the SJ data, the inclusion of the interaction terms GxE and AxE to models ELG and ELA induced a reduction of approximately 10% and 8.7% in the estimated residual variance, respectively. The greatest reduction in residual variance (14.7%) due to interaction effects was observed in the most comprehensive model. A similar trend was observed for all the other traits in SJ and NS data as well. However, the percentage of total variance explained by interaction effect in SJ was slightly higher than that of NS. This indicated that not all the phenotypic variation could be captured fully by main effects alone. A substantial amount of error variance reduction related to the ELG, ELA, and ELGA models could be attributed to genetic-environment interaction, hence the need to factor them into the study.
In the models that included both marker and pedigree effects, the total additive genetic variance was partitioned into fractions due to marker and pedigree information ( Table 4). Both main effect models and interaction effect models were analyzed. The proportion of the additive variance due to the markers and pedigree information was estimated as the ratio of marker (genomic) or pedigree variance to the total additive variance, respectively. In the combined data, most of the total additive genetic variance (average of 60% for the three traits) was attributed to marker effects.
However, a different scenario was observed when data from the different breeding companies were analyzed separately. Averaged across all four traits in the SJ data, 53% of the total additive genetic variance in the main effect model was attributed to the pedigree relationships among the lines. In NS data, however, more than half of the total additive genetic variance for all the traits, with the exception of protein, in the main effect model was attributed to the genomic/marker relationships among the lines. The scenarios observed with the interaction models were similar to the main effect models. This indicated that  both the genomic and pedigree relationship matrices were able to capture significant proportions of the additive genetic variance in the population with the relative importance of genomic and pedigree information varying by trait and breeding company used in the current study.

Heritability
Broad-sense and narrow-sense heritabilities were both estimated from the variance component analysis, with the former based on the sum of the line variance and the pedigree and/or marker variance. Whereas marker and pedigree capture the true additive genetic effect, the lines capture the genetic residual effect (Oakey et al., 2016;Hunt et al., 2018). Broad-sense heritability is needed to adjust prediction ability when computing prediction accuracy, and they can be inferred from the variance components estimates (data not shown). Figure 3 shows the narrow-sense heritability h 2 based on individual line records estimated from the main effect models for all traits for the separate breeding programs. The models compared were those that used marker relationships alone, pedigree relationships alone, and combined marker and pedigree relationships for all the traits.
Overall, h 2 in NS was higher than in SJ for all the traits measured. In SJ, h 2 was highest for SF_abv2.5, followed by SF_abv2.2, test weight, and protein. In NS, h 2 was highest for test weight, followed by SF_abv2.5, protein, and SF_abv2.2. As expected, h 2 traced by pedigree information was higher than that traced by marker information for all the traits in the SJ data. Overall, h 2 was reduced by almost 29, 19, 20, and 25% for protein, test weight, SF_abv2.5, and SF_abv2.2, respectively, when the marker relationship matrix ELG was employed rather than the pedigree relationship matrix (ELA). It was also observed that h 2 traced by combined pedigree and marker information was higher than using either of the matrices on their own.
In NS, however, except for test weight h 2 traced by pedigree information alone was higher than marker information alone   for protein, SF_abv2.5, and SF_abv2.2. In addition, combining pedigree and marker information improved h 2 compared to marker or pedigree information alone for protein and test weight; however, in the seed fractions there was no gain compared to the pedigree information.

Within-Company Predictions Using Different Cross Validations
Prediction results for the different cross validation schemes from the six models using data from the two breeding programs (NS and SJ) separately (within-company predictions) for all traits are presented in Figure 4.

Prediction by Cross-Validation Scheme CV1
In predicting newly developed lines in CV1, pedigree information alone models (ELA) performed better (SJ data) or similar (NS FIGURE 4 | Prediction ability obtained using six models for the three cross-validation schemes in the two breeding programs for all traits. CV1; predicting newly developed lines in known environments, CV0; predicting unknown (future) year, and CV00; predicting newly developed lines in unknown future year prediction abilities for (A) Protein, (B) Test weight, (C) Seed fraction SF_abv2.5, and (D) Seed fraction SF_abv2.2. Error bars depict the standard error from the mean of the 10 fold correlations. data) to marker information alone models (ELG) for all traits. In SJ, prediction ability for ELG in CV1 were 0.27, 0.46, 0.58, and 0.54 for protein, test weight, SF_abv2.5, and SF_abv2.2, respectively, while prediction ability for ELA in CV1 were 0.31, 0.50, 0.60, and 0.58 for protein, test weight, SF_abv2.5, and SF_abv2.2, respectively, (Figure 4). In NS on the other hand prediction ability for ELG in CV1 were 0.44, 0.60, 0.62, and 0.64 for protein, test weight, SF_abv2.5, and SF_abv2.2, respectively, while prediction ability for ELA in CV1 was similar to ELG for all the traits. Combining pedigree and marker information (ELGA) improved prediction ability compared with models using marker information alone, but not compared with models using pedigree information alone for SJ data for all traits. In NS, ELGA models performed similarly to ELG and ELA models for all traits.
Inclusion of interaction effects increased predictions compared to the main effect models for all traits in the CV1 scheme for both SJ and NS. The greater prediction performance from the inclusion of interaction effect with the main effect models was found to be highest for protein in in both SJ and NS. In addition, the gain in prediction performance due to the addition of interaction was better in SJ than NS for most of the traits. For instance, in SJ, GxE improved the prediction ability of the ELG model by approximately 26, 9, 4, and 5% for protein, test weight, SF_abv2.5 and SF_abv2.2, respectively, while AxE improved the ability of the ELA model by 21, 6, 4, and 6% for protein, test weight, SF_abv2.5, and SF_abv2.2, respectively. In NS, inclusion of GxE in the ELG model increased prediction ability by 12, 6, 4, and 6%, while AxE increased the prediction ability of the ELA model by 11, 5, 3, and 5% for protein, test weight, SF_abv2.5, and SF_abv2.2, respectively. Thus, the models that captured the interactions performed better than their corresponding main effect models for all traits.

Prediction by Cross-Validation Schemes CV0 and CV00
The scheme for predicting performance of lines in future year (CV0) and the scheme for predicting the performance of new lines in future year (CV00) showed a different pattern compared to CVI. For all traits in both SJ and NS, marker alone models outperformed pedigree alone models. The average predictive ability for protein content in both companies was 0.19 and 0.21 for marker models in CV0 and CV00, respectively, while with pedigree models; it was 0.07 and 0.03 for CV0 and CV00, respectively. In test weight, the average predictive ability for marker models in both breeding programs was 0.36 (CV0) and 0.34 (CV00) while the average predictive ability for pedigree models was 0.30 (CV0) and 0.18 (CV00). The trend was similar for the other traits.
In NS, though ELGA models performed better than ELA models, there was no advantage in predictive ability compared to ELG models for all the traits under CV0 and CV00 schemes. In SJ on the other hand, ELGA always performed better than ELA but its performance against ELG was inconsistent depending on the trait. In protein, ELGA performed worse than ELG while in test weight, ELGA performed similar to ELG for both CV0 and CV00. In SF_abv2.5 and SF_abv2.2, ELGA performed similar to ELG in CV0 but better in CV00. Interestingly in CV0 and CV00, inclusion of interaction effects did not always improve prediction performance as seen in the case of CV1 (Table 5). In most cases, the addition of interaction effect in CV0 and CV00 caused a decrease in prediction performance. In cases where there was an increase due to interaction, the gains were always lower compared to CV1 with a few exceptions in SF_abv2.5 and SF_abv2.2. For instance, in SJ the percentage change in predictive ability were between approximately -158% (CV00 AxE) to 10% (CV00 GxE) while that in CV1 ranged between approximately 21% (AxE) to 33% (GxE-AxE). In NS, the percentage change in predictive ability due to interaction ranged between approximately 11% (AxE) to 15% (GxE-AxE) for protein content in CV1 and between -6% (CV00 AxE) to 5% (CV00 GxE) for the forward predictions.

Comparing Within-Company and Across-Company Predictions Using CV1
The effect of combining populations to increase training set size in prediction was assessed by comparing predictive ability estimated using within-company (lines from either one breeding program) or across-company (lines from both breeding programs). Training sets that included lines from combined breeding programs increased the entries in that training set by almost twice compared to training sets that included lines from a single breeding company. The validation sets always contained lines from either of the breeding programs. The hypothesis here was that training sets from combined populations would increase prediction performance compared to training sets from single populations since by combining populations, the size of the training set is increased.
Our results were mostly population and trait specific. When the validation is from SJ, then combining populations from NS, and SJ either did not change or slightly increased prediction ability in most of the models (except for the most comprehensive model) for protein and test weight. In SF_abv2.5 and SF_abv2.2, increasing training set size by combining populations from SJ and NS reduced prediction ability in most of models ( Figure 5). When the validation set is from NS, increasing training set size by combining SJ and NS reduced prediction ability in all the traits for all the models. In the most comprehensive model ELGA-GxA-AxE, there was a consistent reduction in prediction ability by combining populations from both breeding programs to predict performance of lines from either of the programs. With the exception of ELGA-GxA-AxE, performance of the other models was similar to the results obtained from the withincompany prediction.

Variance Components and Heritability
This study explored several models to determine the effect of independent and combined marker and pedigree relationships on prediction, as well as the effect of environmental interaction terms. From the results of the variance component analysis, the inclusion of the interaction of markers and environments or pedigree and environments was found to induce a remarkable reduction in unexplained variance for all traits. A similar decline has also been seen in previous studies in wheat (Burgueño et al., 2012;Jarquín et al., 2014;Sukumaran et al., 2017) and cotton (Pérez-Rodríguez et al., 2015). It was also observed that interaction terms contributed to the total within-environment phenotypic variance to differing degrees. This indicated that models with only main effects were not adequate to capture most of the variance in the model. Hence, the introduction of interaction terms meant that an increased proportion of the model variance could be captured, which is likely to increase prediction ability.
Models using both marker and pedigree information increased heritability for most of the traits compared to models with either marker or pedigree information alone. Similar results were also observed in sorghum (Velazco et al., 2019) where a k matrix composed of combined pedigree and marker information increased trait heritabilities. Increased heritability due to combined marker and pedigree relationships indicated that such models have a tendency to capture greater genetic variation than using either marker or pedigree information alone.
The reduction in overall narrow-sense heritabilities when ELG was employed compared with ELA in this study reflects the notion of a pedigree-based relationship matrix not adjusting for Mendelian sampling, thus leading to a likely overestimation of additive genetic variance. The general consensus in the literature is that marker-derived relationship information is evolving to becoming preferred in genomic prediction research to pedigreederived relationship information (Crossa et al., 2010;Albrecht et al., 2011;Burgueño et al., 2012;Auinger et al., 2016;Velazco et al., 2019). This is because, unlike pedigree-derived relationship information, marker-derived relationship information accounts for common ancestry, which is otherwise not considered in the pedigree. Moreover, the marker-derived matrix particularly accounts for departures of realized genetic relationships from their expected values due to Mendelian segregation (Burgueño et al., 2012;Velazco et al., 2019), and thus providing relatively accurate genetic variance components compared with the pedigree-derived matrix (El-Dien et al., 2016).

Model Performance Under Different Cross-Validation Schemes
The trend for model performance in terms of prediction ability was similar for all the traits with the different schemes exhibiting differences in how the models performed. Results from this study showed that model performance was dependent on the crossvalidation scheme used. In general, prediction of NS lines was better than SJ lines for most of the traits. This was driven by the close relationship observed between lines in the matrices for NS, as shown in Figure 2. With a low number of subfamilies within NS compared to SJ, the number of full siblings within each subfamily unit in NS was relatively high compared to SJ, thus contributing to lines from NS being predicted more accurately than SJ. This was to be expected because NS has a more recent barley-breeding nursey with narrower genetic base comprised mostly of Northern Europe varieties, while SJ barley breeding has a much longer history with a broader genetic base consisting of genotypes originating from all over Europe.
In CV1, pedigree information was superior or competitive with marker information. In SJ data, pedigree models performed better than marker models while in NS, pedigree models performed similarly to marker models. This result is consistent with previous observations in wheat and sorghum (Juliana et al., 2017;Hunt et al., 2018;Howard et al., 2019). CV0 and CV00 involve more complex prediction situations that span across different environments (years), and in these situations marker information models performed better than pedigree information models. In these complex situations, the number of related lines (full and half siblings) in the training set reduces when sampling across years, and hence a clear advantage of the realized relationship captured by marker is seen (Juliana et al., 2018). When exceptional pedigree records are kept that go back several generations, pedigree-based models may perform better than marker-based models (Juliana et al., 2017). However, this will only hold if the number of related lines in the training and validation sets is optimal. Thus, when lines are derived from a pedigree that contains diverse siblings, they are more likely to be predicted accurately using marker information than using pedigree information, in which case all full siblings within a single subfamily would have the same predictions (Hunt et al., 2018).
Dense marker coverage is a prerequisite for improved prediction accuracies (Heffner et al., 2009). However, the amount of additive genetic variance that can be explained by the markers not only depends on dense marker coverage, but also on the number of markers in linkage disequilibrium (LD) with causative genes (Jensen et al., 2012). This implies that when markers in LD with causative genes are low, models based solely on markers may perform worse than their corresponding models based only on pedigree. Although it was assumed that the number of markers used in this study provided excellent genome-wide coverage, there is a possibility that these markers insufficiently covered some major regions associated with the traits, thus leading to lower prediction accuracies in some of the models based only on markers. This was evident by most of the total additive genetic variance being explained by pedigree relationships and not marker relationships (Table 4) when the models included both matrices especially in the SJ data.
Combining marker and pedigree information did not show any advantage in prediction ability. Several studies have reported the improvement of prediction by a combination of marker and pedigree relationships with respect to models with marker or pedigree information alone (Crossa et al., 2010;Albrecht et al., 2011;Burgueño et al., 2012;Auinger et al., 2016;Juliana et al., 2017;Velazco et al., 2019). The reported improvement in prediction of combined marker and pedigree information models is based on the assumption that if markers are not in full LD with QTLs, the addition of pedigree relationship information will contribute to capturing associations between causative alleles because of common ancestry identity, which can improve prediction accuracies (Velazco et al., 2019). However, in the presence of high-quality genomic markers and pedigree information, there is a possibility of some redundancy between the regression on markers and that of pedigree, leading to no or only a marginal gain (Albrecht et al., 2011;Juliana et al., 2017). The lack of advantage when both matrices were combined in this study could be attributed to redundancy in the information captured by both  or due to issues of collinearity arising from them since they are not mutually orthogonal (Howard et al., 2019).
Assessing and quantifying genotype-by-environment interaction is an integral part of GS and has been shown to lead to increase in prediction ability . In this study, inclusion of genotype-by-environment interaction terms to the main effect models led to increase in prediction abilities for all the traits in CV1. Protein, which had the most dramatic reduction in within-environment error variance following the addition of interaction terms, had the most dramatic gain from interaction terms with respect to prediction ability in both NS and SJ. This gain in prediction ability was not surprising, given the contribution of the interaction terms to the main effect models in the variance component analysis, as discussed previously. Moreover, in CV1, the most comprehensive model (ELGA-GxE-AxE) had the highest prediction ability across all traits. Bearing in mind the extent of the GxE and AxE variance contribution in each of the traits, this increase was not a surprise.
Thus, it appears that if combining marker and pedigree is favorable, sizable gains in prediction ability can be attained by jointly considering their environmental interactions in barley, as has also been observed in previous studies in wheat (Burgueño et al., 2012;Sukumaran et al., 2017).
Interestingly, accounting for genotype-by-environment interaction did not always improve prediction ability in CV0 and CV00. Similar results have been reported in wheat (Howard et al., 2019) and this could be because CV1 which involves the prediction of newly developed lines in tested environments allows borrowing of information from other environments (years) while CV0 and CV00 which involves the prediction of performance of unobserved environments, a much more complex situation does not.

Increasing Training Set Size by Combining Populations
Several studies have shown that the increasing training set size improves prediction ability (Zhong et al., 2009;Asoro et al., 2011;Zhang et al., 2017). Due to this notion, it was hypothesized that increasing the training set size by combining lines from the different breeding programs would improve prediction ability.
The results indicated that increasing training set size by combining different populations produced no benefit in terms of predicting each population test set. In most cases, it led to either a no change in prediction ability or a decrease. The loss in prediction ability of across-company compared to withincompany ranged from 0.7 to 4.6% on average when the validation set comprised of SJ data and from 3.5 to 4.7% when the validation set comprised NS data. This finding is similar to what has been observed in GS for Fusarium head blight resistance in six-rowed barley (Lorenz et al., 2012) and GS for grain yield in bread wheat (Juliana et al., 2018). Conclusions from both studies suggest that the lack of advantage in combining populations corresponding to different breeding programs to increase the training population could be due to the possibility of marker effects differing across the populations and hence causing substantial QTL-byenvironment interaction. Other reasons could be due to different QTLs segregating between the two populations and different QTL effects between the different populations due to epistasis (Lorenz et al., 2012).
In contrast to these results, others have shown that prediction ability can be increased by combining multiple populations. In animal breeding, genomic selection has been shown to benefit from combining populations when the genetic distance between the breeds is close. In distance related breeds, genomic prediction can be improved if the models used are those that put more focus on genomic markers in strong LD with causative variants (Lund et al., 2014). Similarly, the advantage from combining populations can surge with increasing number of markers, increasing heritability and decreasing divergence time between the populations. With an increased marker density, even very diverged populations are likely to share marker-QTL LD structure that was generated due to mutational events that occurred before the divergence of the population (de Roos et al., 2009). This in effect reduces the negative influence of contrasting QTL-marker effects between the populations and subsequently reduces recombination frequency between markers and QTLs.
The populations used in this study were assumed to be quite related since they come from competing breeding companies in close proximity targeting the agro-climatic zones of Northern Europe. There is high probability of both companies sharing common lines as founders since they have the same consumer target and the possibility of varieties released by one company being used as a source of new genetic material for crosses in the population of the other (Robertsen et al., 2019). However, from our results we can deduce that the populations are not related enough to benefit from this combination. By exploring more sophisticated variable selection methods, the advantage of increasing training population size by combining lines from both breeding programs can be realized, since such models are able to capture tightly linked QTL-marker effects.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
TA-Y, SR, and LJ conceptualized the idea for the manuscript. TA-Y analyzed the data and wrote the manuscript. LJ supported in the statistical analysis. RH and JJ designed the field trials and provided genomic and pedigree data. RH, JJ, and TA-Y collected the phenotypic data. All authors read and approved the final manuscript.