Predictive assessment of single-step BLUP with linear and non-linear similarity RKHS kernels

Abstract


Introduction
Genomic selection (GS)is widely used across plant and livestock species and has been well accepted by genetic improvement companies.GS uses genomic information like single nucleotide polymorphism (SNPs) data, to estimate genomic breeding values and rank selection candidates in a breeding program (Pryce & Haile-Mariam, 2020;VanRaden, 2020).Different statistical approaches and strategies have been used to predict genomic estimated breeding values, GEBV (e.g., Gianola & Rosa, 2015).The most commonly used method based on genomic relationships or similarities (Gianola et al., 2020), is known as genomic best linear unbiased prediction (GBLUP).The method is a modification of traditional pedigree-based best linear unbiased prediction (ABLUP), a standard for predicting breeding values using expected relatedness among individuals derived from pedigree information.GBLUP differs from ABLUP in that the relationship matrix A, is replaced by a genomic relationship matrix (G) that is calculated from genotypic data to capture realized relatedness resulting from the process of Mendelian sampling (Bernardo, 1994;Misztal et al., 2020;VanRaden, 2008).
An important development took place when GBLUP was extended to the "single-step" GBLUP method (ssGBLUP), which allows incorporation of both pedigree-and genomic-derived relationships into a single relationship matrix H (e.g., Misztal et al., 2009).An essential component of the single-step method is that the genomic relationship matrix among genotyped animals is expanded via using pedigree information to form a relationship matrix for all animals, including individuals that were not genotyped.The combined relationship matrix (H) provides a framework for obtaining GEBV of all individuals in the pedigree simultaneously in a single step ( Christensen & Lund, 2010).Early attempts at combining GEBV and breeding values (EBV) were based on blending "direct" genomic values (DGVs) based solely on genomic and phenotype information, with EBVs by using indexes that weighted the two estimates of breeding values in some manner.
Blending DGVs and EBVs was based on the rationale that, if the effect of quantitative trait loci (QTL) was not fully captured by the genomic markers, it could still be captured by polygenic effects (Konstantinov & Hayes, 2010;Pryce & Haile-Mariam, 2020;VanRaden, 2008).In ssGBLUP (Aguilar et al., 2011;Christensen & Lund, 2010), pedigree, phenotypes, and genotypes are used jointly to predict genomic estimated breeding values (GEBVs), for all individuals by using what is essentially an imputation of genomic values using pedigree information that connects genotyped individuals with individuals without genotypes.
The concept of ssGBLUP is operationally attractive because it allows exploiting available computing strategies suited to large-scale BLUP implementations.A number of studies based on either real or simulated data has indicated that ssGBLUP is effective and that predictions can be better than those delivered by DGV or blending methodologies (Howard et al., 2014;Konstantinov & Hayes, 2010;Pérez-Rodríguez et al., 2012).Methods using either SNP effects or genomic relationships were initially based on a multistep approach (VanRaden, 2008), where a regular genetic evaluation by pedigree BLUP was followed by extraction of pseudo-phenotypes for genotyped animals, followed by an evaluation of genotyped animals, and then by a calculation of an index that combined pedigree and genome-based information ( VanRaden, 2008).
On the other hand, several studies suggested that non-parametric methods based on kernels, such as reproducing kernel Hilbert space regression (RKHS) improve predictions of complex traits (Gianola et al., 2006;Gianola & Van Kaam, 2008).In particular, it was conjectured that non-linear Gaussian kernels (GK) could capture complex non-additive gene action (e.g., gene×gene epistatic interactions), as well as nonlinear relations between phenotypes and genotypes.Subsequently, de los Campos et al., (2009), de los Campos et al., (2010), and Pérez-Rodríguez et al., (2012) noted that BLUP or GBLUP are special cases of RKHS.Many studies have suggested that various kernels derived from marker information, could outperform the predictions delivered by the G relationship matrix (González-Camacho et al., 2012;Pérez-Rodríguez et al., 2012), which is a valid kernel for RKHS as well, as noted above.It appears that RKHS can improve prediction accuracy, particularly if there are genotype by environment interaction, epigenetic or metagenomic effects (Cuevas et al., 2016;E Sousa et al., 2017).Cuevas et al.(2019) recently introduced a positive-definite arc-cosine deep kernel (DK) for genomic prediction as an alternative to deep learning (DL) methods, and which retains the theoretical appeal of RKHS of capturing relationships or similarities between individuals.Crossa et al., 2019aCrossa et al., , 2019b, reported that DK achieved a similar or slightly higher prediction accuracy than either the GK kernel or the genomic relationship matrix (G).The tuning parameter "number of layers" required for DK can be found using a maximum marginal likelihood procedure (Cuevas et al., 2019).The number and kind of genotyped individuals are crucial for a successful application of ssGBLUP approach, and these factors impact prediction accuracy in a breeding program (Auinger et al., 2021;Gianola, 2021).For example, a dairy cattle study by Granado-Tajada et al., (2021) using the ssGBLUP approach found that genotyping males and female are beneficial, when these animals possess daughters with lactation records.There was no gain in prediction accuracy when the genetically best (putatively) or extreme individuals were genotyped.They also emphasized the importance of genotyping individuals from several generations.There seems to be little recognition that kernel methods can also be used in single step strategies.
In an attempt to examine their performance in a ssGBLUP setting for genomic prediction, we carried out an experimental comparison using a real chicken data.Our study evaluated the predictive ability under different genotyping strategies of: 1) the Gaussian nonlinear kernel (GK) suggested by Gianola et al., (2006) and 2) an arc-cosine deep kernel (DK) suggested by Cho & Saul (2009), where the kernel evaluation in one of several layers is a function of the kernel values in the previous layer.The methods were compared with the standard ssGBLUP method, with pedigree and genomic relationships, which served as a benchmark.

Data
The dataset used was obtained from Aviagen Ltd (Aviagen Ltd, Newbridge, UK), a major poultry breeding company.The phenotypic measurements considered were body weight at 35 days of age (BW) and hen-house production (HHP, the total number of eggs laid between weeks 28 and 54), with heritability of 0.33 and 0.19, respectively ( Momen et al., 2017).All individuals ( = 5500) were sampled from a broiler chicken line undergoing several generations of selection.Pedigree, genotype and phenotype data were available for all birds.Features and pedigree structure of the dataset are shown in Table 1.
When dealing with polygenic traits and single step genomic prediction model, it is typical in animal breeding genetic evaluation that there will be a large pedigree, involving both genotyped and nongenotyped individuals.Use of all available information is desirable to mitigate biases due to selection (Im et al., 1989).Poultry breeding programs maintain complete and deep pedigrees of all birds and employ BLUP-for predicting breeding values of all selection candidates.However, the number of genotyped individuals typically smaller than the number of individuals in the analysis.
< Table 1 about here> Phenotype correction and genotype quality control Before assessing the predictive performance under different design scenarios of the various models, the raw phenotypes were corrected for plausible (treated as fixed) environmental effects, to remove known nuisance non-genetic sources of variation.
All birds were genotyped using a 50K SNP panel from ThermoFisher.Quality control consisted of eliminating SNPs with a minor allele frequency lower than 1% (MAF<0.01)and a call frequency lower than 0.95.A total of 42,780 SNPs remained for downstream analysis after the quality control.

Statistical Analysis
We considered a single trait -single step BLUP model, that included both marker and pedigree information simultaneously for computing the genetic evaluations.Following Legarra et al., (2009) and Aguilar et al., (2010), the model and related variance-covariance matrices were : where  is the vector of corrected phenotypes); 1 is a vector of ones, and μ is the overall mean.Z is the incidence matrix that related observations to random genetic additive effects.Term u is the vector of random genetic additive effects, is assumed to follow the multivariate normal distribution (,  !" ), where  !" is the variance of additive genetic effects; e is the vector of random residual effects, following the normal distribution of (,  # " ), in which  # " is the residual variance.Let  = ( $ % , ′ " )' where the partitions pertain to non-genotyped and genotyped individuals, respectively.The matrix H was defined as: In the above expression, A11, A12, A21 and A22 are sub-matrices of A (the pedigree-based relationship matrix))Here, can be any  ×  positive-definite kernel matrix which reflects the covariance structure (i.e., conveying molecular similarity) between the genotyped individuals.The kernel matrix K is built from marker information, using various operations on marker codes.

Designing the proportion of genotyped individuals
In our dataset there were genotypes for all individuals.To mimic the single step BLUP setting comprising genotyped and non-genotyped subsets, we designed three genotyping scenarios by masking a varying portion of the entire marker genotypes of individual birds with pedigrees, to construct the H matrix.In each of the three scenarios 20%, 40%, 60% or 80% of individuals had genotypes.For example, in the 20% setting, all birds had pedigrees but only 20% were presented to the model with marker information.
The first scenario was called "youngest individuals genotyped" (YS).Here, kept the genotypic information on subsets of 20% ( = 1100), 40% ( = 2200), 60% ( = 3300) and 80% ( = 4400) according to age of the bird.For instance, in the 40% setting, the 40% youngest birds in the pedigree were presented with genotypes.These individuals (in the 40%) had all phenotypic, genotypic and pedigree information, and the rest of the birds had only pedigree and phenotypic information, so their genotypes were masked.
In the second scenario (PA), individuals with genotype information were selected based on average phenotype of its parents.The parental average for each bird was calculated  &'()#*+ = 0.5( ,-'# +  ./0), where y is the adjusted phenotype.Individuals with missing information for both parents were discarded; if there was information on only one of the parents, the evaluation was the adjusted record of the single parent.We selected 20%, 40%, 60% and 80% of the top averages as individuals possessing genotypic information.For example, in the 80% genotyping setting, only 20% of the birds had pedigree data and phenotypes, whereas 80% had marker, pedigree and phenotypic data.
Finally, for the third scenario (RS), we randomly selected subsets of 20, 40, 60, and 80% of the individuals from all genotyped animals in the dataset sets.In contrast to the two previous scenarios, no consideration of age or performance was made in this scenario.

Kernel methods
We constructed three different similarity kernels based on the additive encoding of marker effects (K); these kernels were then used in the RKHS regression model (de los Campos et al., 2009).The first kernel was the genomic relationship matrix suggested by VanRaden (2008), typically used in ssGBLUP and applied to various species: where, M is a n×m centered genotype incidence matrix for individuals (: 1, 2, … , ) of SNP additive codes (j=1, 2,…, m; m= number of markers).SNP genotype codes were  -2 ∈ {0 = ; 1 =   ; 2 = } and p 1 is the allelic frequency of the minor allele at the j-th SNP.
The G matrix was used to construct H in a single step BLUP that was used as benchmark for comparisons.

Gaussian kernel
The nonlinear Gaussian Kernel (GK) method (e.g., Gianola et al. 2006;Gianola and van Kaam 2008) was the second type of kernel used.The Gaussian kernel has the form: where ‖ 4 - 4 $ ‖ is the Euclidean distance between the vectors of SNP markers of individuals i and i % normalized to range from 0 to 1, relative to the median of the pairwise distance Q, a scalar variable; h > 0, is a bandwidth parameter (a regularization variable) that controls the similarity between individuals or rate of decay of GK 44 $ Euclidean distance (increase or decrease).GK, is one of the most widely used kernel functions in genome-enabled prediction, and selection of the bandwidth is critical.Here, we used an approach called "kernel averaging" or "multiple kernel learning," as proposed in de los Campos et al., (2010).We defined a grid of seven values: h= (0.2, 0.4, 0.8, 1, 1.5, 3, 5) and using the formula above, we computed seven distinct GKs, named GK0.2,GK0.4,…, GK5, related to each specific value of h.Then, the seven kernels were "averaged" to build a final kernel as:  = 5 %& '.) ) , …,  ;< + " parameters are variance component estimates captured by the kernels  7." ,  7.9 ,…,  : , respectively, and  ] ;< " is the sum of these seven variances.We assumed that the ratios of variances reflect the relative contributions of the kernels to the marked genetic variation in the population (Supplementary Excel spreadsheet).The resultant  kernel matrix was used to construct  to be used in a single step GBLUP.

Deep Kernel
The arc-cosine kernel, referred to as Deep Kernel (DK), was the third similarity matrix employed (1 + ())be the "Heavyside" step function taking the value zero for negative arguments and one for positive arguments.We defined the t-th order of arc-cosine kernel function by integral representation: where, w is the weight corresponding to the parameters of the model.For non-negative integer values of , Cho (2012), showed that the the angle  between the  4 and  4 $ input vectors is: Here, ⋅ stands for the inner product and ǁMiǁ is the norm of individual's i genotypes.The kernel resulting from the above operation is a symmetric semi-definite positive matrix (Cuevas et al. 2019a).For a single layer in an artificial neural network (ANN) layout, let: Where, π is the pi constant and r -,- takes its maximum value at  = 0, and decays monotonically to zero at  = p, for all values of t.
When  = 0, the arc-cosine kernel maps inputs M, to a unit hypersphere in feature space with  7 (, ) = 1;when  = 1, the arc-cosine kernel preserves the norm of inputs as  $ (, ) = ‖‖ " .Finally, for all  > 1, the kernel is A K (, ) ~ ‖‖ "= .A potential advantage of DK is the ability of capturing non-additive relationships between individuals, an unexplored concept in quantitative genetics theory.Cho and Saul (2009) and Cuevas et al. (2019) present a recursive relationship approach for shaping a basic DK, into a final DK-emulating ANN hidden layer (), by repeating  times the operation Thus, values of  (NF$) , at level (layer)  + 1 are computed from the previous layer  (N) .
Computing a bandwidth is not necessary, contrary to GK, and the additional computational effort required depends on the number of discrete layers.We selected the number of layers (l), using the maximum likelihood method in (Cuevas et al., 2019).

Prediction ability by cross-validation (CV)
Genome ), where r stands for rank correlation, to normalize the distribution of correlation estimates.We also performed a test for empirical prediction bias done by regressing phenotypes on predicted genetic values; if the slope of the regression differs from 1, this would suggest "bias".
All models were fitted with the "emmreml" function from the EMMREML R package (Akdemir and Godfrey 2015).

Predictive performance for body weight at 35 days of age (BW)
Figure 1 shows the boxplot of the predictive rank correlations, PMSE and bias (assessment (slopes) values for the three genotyping scenarios (PA, RS and YS) for body weight (BW) over the 200 replicates.In all genotyping scenarios, recall that we first selected 20% of the youngest genotyped chickens as animals with genotypic information in the H matrix.Then, we allowed to have genotypes to 40% , 60% and 80% of the birds in the sample (percentage values on the x-axis).The predictive performance of different H matrices is indicated by blue, red and green colors, for G, AK and DK, respectively.Predictive rank correlations increased as the proportion of birds with genotypes increased from 20 to 80 %.This was observed for all three H matrices, across all genotyping scenarios (PA, RS and YS).For example, in the first column of Figure 1 (PA), for the scenario with 20% genotyped birds, the mean predictive rank correlations (standard deviation) were 0.25 (0.02), 0.25 (0.02) and 0.26 (0.02), for HG, HAK and HDK, respectively, and increased to 0.31 (0.02), 0.31 (0.02) and 0.33 (0.02) when 80% of birds were genotyped.Under the 80% setting, there was a mild advantage of HDK over HG and HAK in the single step BLUP models.In YS, birds selected for genotyping according to their age in the pedigree, the most closely related animals originated from recent generations.YS is representative of a selection scenario where genotyping and phenotyping of youngest progenies is favored.Here, there was much overlap between the predictive distributions generated by different kernels, with slight advantage for HDK.
A similar pattern of mild differences between kernels was observed for predictive mean squared error (PMSE) for all genotyping scenarios.As the fraction of genotyped individuals relative to the total increased, PMSE decreased; the lowest PMSEs were obtained with 80% genotyped individuals with HDK.As depicted by the bottom plots of Figure 1, the slopes could not be considered different from 1, so all predictions could be claims empirically "unbiased".
The panels in the middle column in Figure 1 compare HAK and HDK versus HG when individuals were randomly genotyped (RS).Under RS, all three kernels for single step BLUP, had poorer prediction ability when compared to the YS and PA scenarios.Ass before, the lowest and highest prediction rank correlations were obtained with 20% and 80%, genotyping, respectively.The ranges of predictive correlations under RS were 0.10 (0.02), 0.11 (0.02), and 0.10 (0.02) with 20 % genotyped individuals, and increased to 0.15 (0.02), 0.18 (0.02) and 0.17 (0.02) with 80 % genotyped birds, for HG, HAK and HDK, respectively.There was a hint of a superiority of HAK and HDK, over HG but it did not translate into lower MSE.The leftmost column of Figure 1 shows the predictive performance of the H matrices when birds were selected for genotyping based on the phenotypic parent average (PA).For this scenario, the lowest predictive rank correlations were again obtained when only 20 % of the birds were genotyped, with values 0.24 (0.02), 0.25 (0.02) and 0.25 (0.02) for HG, HAK and HDK, respectively; the largest values were obtained with 80 % of individuals genotyped.In the PA scenario, the HDK, was slightly better than HG and HAK, except when only 20 % of the birds were genotyped.In short, for PA, HDK and HAK were slightly better than HG.Predictions were empirically "unbiased" in YS since the slopes of the regressions did not differ from 1. Overall, predictions were better in the YS and PA scenarios and worst in RS in terms of all metrics considered.
In a nutshell, results body weight (BW) indicated that single step BLUP predictions may be improved in some cases by using non-linear similarity matrices for the H matrix, without detectable adverse effects.This result held mostly when predictions derived from a large proportion of individuals with genomic data, in addition to pedigree and phenotypic information.
The non-parametric kernels have potential to capture additive and non-additive gene actions (Morota & Gianola, 2014), and this property is expected to be conveyed to some extent to the H matrix.In general markers exploit similarity in state, and may capture non-additive gene action (if appropriately encoded) and linkage disequilibrium, whereas A informs about similarity by descent ( Momen et al., 2017), so there would be complementarity between genomic and pedigree data.
The additive encoding of markers and the standard genomic relationship matrix are supplementary to the information from A. Our findings suggest that H matrices employing nonlinear kernels may be useful for attaining a higher accuracy of predictions, when non-additive genetic variance is present without a deterioration in the capture of additive effects, at least in the sense of prediction.genotyping rate was the highest (80 %).HAK was the worst performer under all genotyping rates in RS.Rank correlations ranged from 0.08 (0.02), 0.06 (0.02), and 0.07 (0.02) for 20 % genotyping rate to 0.20 (0.02), 0.19 (0.02), and 0.21 (0.02) for the 80 % rate, for HG, HAK and HDK, respectively.As for BW, RS delivered the lowest predictive ability for HHP.This is in agreement with the view that genomic prediction of more closely related genotyped individuals would be better than of a randomly sampled set of individuals (Pszczola et al., 2011).

<<
In the PA scenario, HAK performed better than HG, and HDK at all genotyping rates.HG had the lowest performances at all genotyping rates.A negligible difference was observed in the predictive rank correlations at 80 % genotyping.The predictive rank correlations were 0.18 (0.02), 0.20 (0.01) and 0.19 (0.02) with 20 % genotyping rate and 0.23 (0.02), 0.24 (0.02), and 0.24 (0.02) with 80 % genotyping rate, respectively for HG, HAK and HDK.Predictive mean squared error (PMSE) displayed the same pattern as predictive rank correlations but differences were minor.No evidence of empirical "bias" was detected.
In summary, predictive accuracies of single step genomic prediction based on non-linear similarity matrices were slightly better, seldom worse, than those based on the traditional single step GBLUP (HG) for BW and HHP.Sometimes imperfect LD can lead to apparent epistasis.A recent study by Schrauf et al., (2020) provided evidence that at a higher marker density the superiority of nonlinear over the standard additive kernel may dissipate if such phantom epistasis exist.In practice, however, it is almost impossible to claim that LD between markers and quantitative trait loci and markers is "perfect" or "imperfect", as the true LD cannot be observed.
When a larger proportion of birds with genotypic information was present in the reference population, the gain of nonlinear kernels was larger, and was larger when genotyping was based on PA, especially when we used HAK.

<< Figure 1 About Here >> Discussion
We investigated two non-linear similarity kinship matrices, the averaged Gaussian kernel (AK) and arc-cosine kernel referred as the deep kernel (DK), when constructing the H matrix in the single-step BLUP methodology.The predictive ability of these kernels was compared with results from a standard genomic relationship-based H matrix (Christensen & Lund, 2010;Legarra et al., 2009).We employed four different genotyping rates, with 20, 40, 60 and 80 % of birds in a sample of fully pedigreed commercial broiler chickens genotyped and, in all cases, the birds had complete phenotype and pedigree information.A training-testing layout was used in every instance and was repeated 200 times by random reconstruction in each scenario of genotyping rate and genotyping strategy (random selection, and based on either age or parental average).
The predictive ability of the models was assessed by comparing the distributions of predictive rank correlations, predictive mean squared errors and prediction bias statistics, and the target traits were body weight at 35 days of age (BW) and hen house egg production (HHP).The latter is the total number of eggs per hen laid between weeks 28 and 54 of age.The two traits have a moderate genomic heritability between 0.19 and 0.36 (Momen et al., 2017), with a negative genetic correlation between them ( )(OP,QQR) = −0.2,Momen et al., 2017).Since all birds had genotype and pedigree information, in order to mimic the setting of single-step BLUP methodology, the genotypes of 80, 60, 40 and 20% of the birds were masked to created varying genotyping rates.
We designed three strategies to decide which genotypes would be masked, to create marker-based kinship matrices and, subsequently, the corresponding H matrices.In the first strategy (YS), we sorted individuals from oldest to youngest, and selected 20, 40, 60 and 80% of the youngest birds.
On a second strategy, we sorted birds according to their parent's phenotypic average (PA) and kept the genotypes of individuals with the highest PA values, with the same as proportions considered for two other strategies (i.e., 20, 40, 60 and 80 %), and finally, we randomly masked genotypes of 20, 40, 60 and 80% of birds (RS).
In all three selection strategies (YS, RS and PA), higher predictive accuracies were obtained for BW than for HHP, as expected based on heritability values of these two traits (Momen et al., 2017).
Predictive accuracy was clearly influenced by the proportion of genotyped birds, with a higher proportion of genotyped birds resulting in a higher prediction accuracy of genomic values.Our results agreed with findings of Boligon et al., (2012) and Chu et al., (2018), who found that selective genotyping improved the accuracy of GEBV, and that animals with the best performance were the most informative for prediction.Selective genotyping is feasible in broilers because important traits such as body weight and feed efficiency can be measured before sexual maturity.
A simulation study by Ehsani et al., (2010) reported that selective genotyping of the best animals based on phenotypic values provided weaker predictions of breeding values of animals in the next generation relative to random sampling, which does not agree with our real data results.In addition, Ehsani et al., (2010), did not find relevant differences between genotyping individuals with high phenotypic values versus individuals with low phenotypic values in the reference population.We found that selective genotyping according to parent average (PA) may deliver a higher prediction accuracy.Using a simulation, Jiménez-Montero et al., (2012) concluded that the predictive accuracy of GEBV depends not only on the number of animals genotyped but also on the selective genotyping strategy as well.
Many efforts have been conducted to enrich BLUP by using alternative kinship-based prediction methods.The significant work of Misztal et al., (2009), well known as ssGBLUP, provided the ability to evaluate genotyped and non-genotyped individuals simultaneously.The methodology has been mostly used for large field data sets, e.g., cattle, pigs and chickens, has led to a higher accuracy, and is simpler than multistep genomic selection methods (Aguilar et al., 2011;Christensen et al., 2012;Simeone et al., 2012).On the other hand, different linear and non-linear marker-based similarity matrices have been developed and implemented by researchers to quantify resemblance between individuals.A commonly used kernel is the Gaussian kernel (GK) based on molecular markers (Gianola et al., 2006;Gianola & Van Kaam, 2008) and recently, Cuevas et al., (2019), introduced the arc-cosine kernel function for genome-enabled prediction.Except for the genomic relationship matrix (a special form of reproducing kernel), none of these kernels have been tested in the context of ssGBLUP.We designed the study to investigate the impact of two well used kernels, AK and DK, in genomic prediction, in a comparison with G, and in the context of ssGBLUP prediction.Our results suggest that, for some of the scenarios tested, the predictive ability of the single step approached can be enhanced somewhat by using an H matrix based on AK and DK, as opposed to G. We found that when genotyping rate increased as part of the selection strategy, the predictive ability of the single step models increased, with the alternative kernels producing in some cases better results than G.Because these kernel methods can capture complex gene actions, as well as nonlinear relationships between phenotype and genotype (Gianola et al., 2006;Gianola & Van Kaam, 2008), the extended H matrix may be useful in predictive problems when dominance and epistasis underlie gene action.Our results also suggest that even when the genotyping rate was small, the prediction accuracy using AK and DK was nearly similar to that of the G, but these two kernels displayed advantages over G at highest genotyping rates.Cuevas et al., (2019), stated that, DK is computationally easier, since no tuning parameter is required, while performing similarly or slightly better than the common kernels.In our study, we used an average of Gaussian kernels, producing a slightly better performance than DK in some cases.A difficulty with the AK approach is that the weights assigned to each of the kernel depend on variance components derived from a multi-kernel fitting exercise.Since the kernels are not mutually orthogonal, the weights placed to the individual kernels may not produce the best possible average kernel.This is a subject for further study.
We used a somewhat large data set representative of poultry breeding studies with genotyped and nongenotypes individuals evaluated together.We found that the AK and DK kernels were slightly better than G, when genotyping rate in the single step strategy was large.Using a high-density SNP panel would be expected to deliver better predictions and perhaps, as in Schrauf et al., (2020), the suggested superiority of nonlinear kernels might be lost as marker density increases, provided that "phantom epistasis") is an illusion created by the LD picture captured by low density panels.The preceding may or may not hold in practice, but our, improvements in prediction accuracies should be taken with caution as evidence of non-additive effects.Through development of new genotyping platforms, the cost of genotyping has steadily decreased, and genotyping a large proportion of individuals will be even more feasible.Whereas computation with large-scale genomic data still remain a challenge, kernel based methods are less involved than marker-based regression approaches.
In conclusion, we studied, seemingly for the first time in the literature, non-linear kernels as an alternative to G in the context of a single step genomic evaluation strategy.The results suggested that this type of kernel may enhance prediction models by capturing additive and non-additive genetic variability, when present.Future research should examine these kernels for traits known to be strongly affected by epistasis or by genotype ´ environment interaction.

Tables
to create an H matrix.The DK structure was introduced inCuevas et al., (2019) and used byCrossa et al., (2019b) for genomic prediction in a multi-environment model.The method is based on Neal (2012), in the context of Bayesian inference for deep artificial neural networks (ANN).An arccosine kernel is used to measure the similarity between two genotyped individuals by considering the angle between two vectors of their SNP markers  -,  -$ .Let Θ() = $ " -enabled prediction accuracy of the various models across the three scenarios was assessed by designing a replicated partitioned training -testing (TRN-TST) layout.Here, training and testing sets in a random partition are completely disjoint.In total, we used 200 TRN-TST replicates, with 60% of the whole data set assigned to TRN and the remaining 40% assigned to TST set in each run.TRN-TST sets were randomly recreated in each replication.The training set was used to fit the models and the testing set to measure the predictive performance of the competingmodels.For each TRN-TST scenario, two metrics were computed: (i) rank correlation between observed phenotypic values and predicted genomic values, and (ii) mean-squared error of prediction (PMSE), i.e., the average squared difference between predicted genomic breeding values and the actual phenotypes.We used Fisher's z-transformation ( % = 0.5[(1 + ) -(1 − )]

Figure 1
About Here >>Predictive performance for hen-house egg production (HHP) Figure2displays the boxplot of rank correlations, PMSE and slope values ("bias" assessment) obtained from the different H matrices over 200 replicates of the TRN-TST layout for hen-house egg production (HHP).Results for YS (right-most column in Figure2), shows a slightly better performance of HDK over HG and HAK when 20 %, 40 % and 60 % of genotyped individuals were used, HDK and HAK kernels had a similar performance for 80 % genotyping rates and HG was slightly worst in this case.The rank correlations for HG, HAK and HDK ranged, respectively, from 0.19 (0.02), 0.19 (0.02), and 0.20 (0.02) for 20 % genotyping rate to 0.21 (0.02), 0.23 (0.02), and 0.23 (0.02) for 80 %.Under the RS scenario, HDK was slightly better than HG and HAK, when

Figure legends: Figure 1 .
Figure legends: Figure 1.Boxplot of the Fisher's z-transformed predictive rank correlations, predictive mean squares errors (PMSE), and predictive bias between phenotypes and predicted breeding values, using H matrices based on the VanRaden's genomic relationship matrix (G), Averaged Gaussian kernel (AK) and Deep kernel (DK) for body weight (BW) in the single step GBLUP model (ssGBLUP).Genotyping scenarios (bottom to top) were: 20, 40, 60 and 80% of birds with genotypes: youngest (YS), at random (RS) and best parent average (PA).Distributions are based on 200 training-testing sets by assigning 60 % and 40 % of birds to training and testing, respectively.Green, red and yellow colors denote values for G, AK or DK relationship matrices, respectively.Figure 2. Boxplot of the Fisher's z-transformed predictive rank correlations, predictive mean squares errors (PMSE), and predictive bias, using H matrices based on the VanRaden's genomic relationship matrix (G), Averaged Gaussian kernel (AK) and Deep kernel (DK) for hen-house production (HHP) in the single step GBLUP model (ssGBLUP).Genotyping scenarios (bottom to top) were: 20, 40, 60 and 80% of birds with genotypes: youngest (YS), at random (RS) and best parent average (PA).Distributions are based on 200 training-testing sets by assigning 60 % and 40 % of birds to training and testing, respectively.Green, red and yellow colors denote values for G, AK or DK relationship matrices, respectively.

Figure 2 .
Figure legends: Figure 1.Boxplot of the Fisher's z-transformed predictive rank correlations, predictive mean squares errors (PMSE), and predictive bias between phenotypes and predicted breeding values, using H matrices based on the VanRaden's genomic relationship matrix (G), Averaged Gaussian kernel (AK) and Deep kernel (DK) for body weight (BW) in the single step GBLUP model (ssGBLUP).Genotyping scenarios (bottom to top) were: 20, 40, 60 and 80% of birds with genotypes: youngest (YS), at random (RS) and best parent average (PA).Distributions are based on 200 training-testing sets by assigning 60 % and 40 % of birds to training and testing, respectively.Green, red and yellow colors denote values for G, AK or DK relationship matrices, respectively.Figure 2. Boxplot of the Fisher's z-transformed predictive rank correlations, predictive mean squares errors (PMSE), and predictive bias, using H matrices based on the VanRaden's genomic relationship matrix (G), Averaged Gaussian kernel (AK) and Deep kernel (DK) for hen-house production (HHP) in the single step GBLUP model (ssGBLUP).Genotyping scenarios (bottom to top) were: 20, 40, 60 and 80% of birds with genotypes: youngest (YS), at random (RS) and best parent average (PA).Distributions are based on 200 training-testing sets by assigning 60 % and 40 % of birds to training and testing, respectively.Green, red and yellow colors denote values for G, AK or DK relationship matrices, respectively.

Table 1
Pedigree information and features of the chicken data used Is a path in the pedigree of an individual, which connects the individuals to its farthest ancestor.