# Deep learning and genomic best linear unbiased prediction integration: An approach to identify potential nonlinear genetic relationships between traits

F. Shokor, P. Croiseau, H. Gangloff, R. Saintilan, T. Tribout, T. Mary-Huard, B.C.D. Cuyabano

<figure><img src="image_1.png" /></figure>

## To cite this version:

F. Shokor, P. Croiseau, H. Gangloff, R. Saintilan, T. Tribout, et al.. Deep learning and genomic best linear unbiased prediction integration: An approach to identify potential nonlinear genetic relationships between traits. Journal of Dairy Science, 2025, 108 (6), pp.6174-6189. 10.3168/jds.2024-26057. hal-05218863

HAL Id: hal-05218863

https://hal.inrae.fr/hal-05218863v1

Submitted on 22 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Deep learning and genomic best linear unbiased prediction integration: An approach to identify potential nonlinear genetic relationships between traits

F. Shokor,<fnref n="1" />,<fnref n="2" />* P. Croiseau,<fnref n="2" /> H. Gangloff,<fnref n="3" /> R. Saintilan,<fnref n="1" />,<fnref n="2" /> T. Tribout,<fnref n="2" /> T. Mary-Huard,<fnref n="3" />,<fnref n="4" />
and B. C. D. Cuyabano<fnref n="2" />

## ABSTRACT

Genomic prediction (GP) aims to predict the breeding values of multiple complex traits, usually assumed to be multivariate normally distributed by the largely used statistical methods, thus imposing linear genetic relationships between traits. Although these methods are valuable for GP, they do not account for potential nonlinear genetic relationships between traits in scenarios. For individual traits, this oversight may minimally affect prediction accuracy, but it can limit genetic progress when selection involves multiple traits. Deep learning (DL) offers a promising alternative for capturing nonlinear genetic relationships due to its ability to identify complex patterns without prior assumptions about the data structure. We proposed a novel hybrid model that combines both DL and GBLUP (DLGBLUP), which uses the output of the traditional GBLUP, and enhances its predicted genetic values (PGV) by accounting for nonlinear genetic relationships between traits using DL. We simulated data with linear and nonlinear genetic relationships between traits in order to verify whether DLGBLUP was able to identify nonlinearity when present and avoid inducing it when absent. We found that DLGBLUP consistently provided more accurate PGV for traits simulated with strong nonlinear genetic relationships, accurately identifying these relationships. Over 7 generations of selection, a greater genetic progress was achieved with PGV that accounted for nonlinear relationships (DLGBLUP), compared with GBLUP. When applied to a real dataset from the French Holstein dairy cattle population, DLGBLUP detected nonlinear genetic relationships between pairs of traits, such as conception rate and protein content, and SCC and fat yield, although, no significant increase in

prediction accuracy was observed. The integration of DL into GP enabled the modeling of nonlinear genetic relationships between traits, a possibility not previously discussed, given the linear nature of GBLUP. The detection of nonlinear genetic relationships between traits in the French Holstein population when using DLGBLUP indicates the presence of such relationships in real breeding data, suggesting that it may be relevant to further explore nonlinear relationships. This possibility of nonlinear genetic relationships between traits offers a different perspective into multitrait evaluations, with potential to further improve selection strategies in commercial livestock breeding programs. This is particularly relevant when integrating new traits into multitrait evaluations or incorporating new subpopulations, which may introduce different forms of nonlinearity. Finally, it is shown that DL can be used as a complement to the statistical methods deployed in routine genetic evaluations, rather than as an alternative, by enhancing their performance.

**Key words:** machine learning, GBLUP, multitrait models, genetic evaluation, genetic relationship

## INTRODUCTION

Genomic prediction (GP; Meuwissen et al., 2001) in genetic evaluations uses DNA marker information, most commonly SNP data, to predict the genetic value (GV) of an individual for complex traits. For each trait, an individual has a single GV, representing the overall genetic contribution to their phenotype. The GV includes both additive and nonadditive genetic effects, but when focusing solely on the additive genetic component that is transmitted to the next generation, it is referred to as the breeding value (BV). In both livestock and plant studies, BV are commonly used for selection purposes. Traditionally, statistical methods applied to GP rely on linear mixed models (VanRaden, 2008), using SNP genotypes, or both joined into the single-step method (Misztal et al.,

<footnote n="1">¹Eliance, 75012 Paris, France</footnote>

<footnote n="2">²Université Paris-Saclay, INRAE, AgroParisTech, GABI, 78350 Jouy-en-Josas, France</footnote>

<footnote n="3">³Université Paris-Saclay, INRAE, AgroParisTech, UMR MIA Paris-Saclay, 91120 Palaiseau, France</footnote>

<footnote n="4">⁴Université Paris-Saclay, INRAE, CNRS, AgroParisTech, GQE–Le Moulon, 91190 Gif-sur-Yvette, France</footnote>
2009). The predicted genetic values (PGV) are obtained under the assumption of normally distributed effects, either through the BLUP or Bayesian approaches, with their various prior assumptions and alphabets (Meuwissen et al., 2001; Erbe et al., 2012). Combined with the current availability of genomic information, these statistical methods are very powerful and have undoubtedly revolutionized genetic evaluations. Due to the assumption of normality of the data, when extended to a multitrait scenario, the relationship between traits is assumed to be linear (Henderson and Quaas, 1976; Calus and Veerkamp, 2011; Montesinos-López et al., 2016). If traits present nonlinear genetic relationships, the prediction accuracy of the current statistical methods may not achieve its full potential (Sölkner and Fuerst-Waltl, 1996; Fuerst-Waltl et al., 1997), despite the efforts of improvement by providing highly informative data, such as denser SNP chips, or using a single-step strategy that combines the pedigree to the genomic data into a blended model (Misztal et al., 2009).

The restriction to linearity, in addition to the continuously increasing amount of recorded and genotyped animals, which posed computational constraints for genetic evaluations with the classical statistical models, were among the reasons that pushed geneticists to explore the use of deep learning (DL; LeCun et al., 2015) for genetic evaluations and genomic predictions (Montesinos-López et al., 2021). Deep learning approaches (LeCun et al., 2015) use artificial neural networks, which have the capacity to learn complex patterns and features from very large datasets and map input to output data. Depending on the prediction task to be achieved, different DL architectures have been considered (see e.g., Vaswani et al., 2017; Shinde and Shah, 2018; Alzubaidi et al., 2021; Jubair et al., 2021; Choudhary et al., 2022).

Most studies that use DL for GP have been motivated by the potential of accounting for nonadditive genetic effects (epistasis and dominance) with the objective of increasing prediction accuracy by accounting for such effects (Bellot et al., 2018; Ma et al., 2018; Pérez-Enciso and Zingaretti, 2019; Abdollahi-Arpanahi et al., 2020; Waldmann et al., 2020; Zingaretti et al., 2020; Pook et al., 2020; Costa et al., 2022; Nazzicari and Biscarini, 2022; Wang et al., 2023; Chen et al., 2024). More recently, Lee et al. (2023) proposed deepGBLUP, a combination of DL and GBLUP into a single model. Given input SNP data, the DL networks of deepGBLUP extract the effects of adjacent SNPs using locally connected layers to estimate a GV through fully connected layers. In parallel, deep-GBLUP fits a GBLUP model accounting for additive, dominance and epistasis effects (through their respective genetic relationship matrices), finally summing all estimated genetic effects (DL, additive, dominance, and epistasis). Applied to a real dataset of Korean native

cattle, and to simulated data, the results of the proposed deepGBLUP model outperformed those of the traditional GBLUP and Bayesian methods in different single-trait prediction scenarios.

In livestock breeding, the goal is to improve multiple commercially important traits simultaneously to achieve overall genetic progress. When working with multiple traits, genetic relationships between traits (van Rheenen et al., 2019) must be considered because selection for one trait will affect the other correlated traits (Hill, 2013). It is thus relevant to account for the influence of genetic relationships on the different stages of a genetic evaluation, from the estimation of variance components to the prediction of GV, in order to jointly improve multiple traits of interest. For example, it has been shown that it is more advantageous to perform a multitrait model for genetic evaluation that explicitly considers the genetic relationships between traits to improve their predictive ability (Montesinos-López et al., 2018; Pérez-Rodríguez and de los Campos, 2022).

Although nonlinear relationships are not straightforward to handle by statistical methods, such a hypothesis is suitable for the implementation of DL methods. Few studies have explored the possibility of nonlinear genetic relationships between traits using statistical methods like polynomial regression in animals (Sölkner and Fuerst-Waltl, 1996; Fuerst-Waltl et al., 1998; Castañeda-Bustos et al., 2017) and humans (Akingbuwa and Nivard, 2023). Even so, such studies focused merely on how to identify these nonlinear genetic relationships, rather than on how to account for them in GP.

The aim of this study was to explore the impact of nonlinear genetic relationships between traits on the performance of GBLUP and DL methods. The specific objectives of our work were to (1) develop a method that uses DL to model these potential nonlinear relationships; (2) compare the accuracy of PGV obtained with GBLUP and with DL methods tailored to identify nonlinear relationships, for traits that present such intertrait nonlinear relationships; (3) evaluate how the presence of these relationships can affect the genetic progress over generations by comparing scenarios where they are ignored (GBLUP) versus where they are accounted for (DL methods); and (4) verify if nonlinear genetic relationships can be identified in a real dataset of the French Holstein dairy cattle population. The GBLUP predictions were taken as the baseline result for comparison, and we proposed a hybrid DL model that combines both DL and GBLUP (DLGBLUP) to leverage their complementary strengths. Although DL can capture nonlinear genetic relationships, which GBLUP cannot, GBLUP provides more accurate predictions of additive genetic values. By integrating both approaches, DLGBLUP aims to enhance prediction accuracy while reducing the limitations of each method.
The DLGBLUP method consists of 2 steps: first, the GV are predicted using a multitrait GBLUP based on the genomic data, and then GV are re-predicted using a DL model that captures potential nonlinear genetic relationships between traits. A full DL model was also studied, in which DL replaced GBLUP for the initial trait prediction. Unlike GBLUP, this approach does not assume only linear genetic relationships. This full DL model was implemented to assess whether replacing GBLUP, which is restricted to linear relationships, with a model capable of capturing nonlinearity could improve prediction accuracy.

# MATERIALS AND METHODS

## Simulated Data

The complete simulated genomic data consisted of 10,000 SNPs distributed across 29 chromosomes, with an average linkage disequilibrium (LD) pattern resembling that of a cattle population (de Roos et al., 2008). From the 10,000 simulated SNPs, 512 were assigned as QTL to be shared between all traits. The noncentered SNPs were coded as 0, 1, and 2, referring to homozygous for reference alleles, heterozygous, and homozygous for alternate alleles, respectively.

A quantitative reference trait was simulated using the 512 QTL for a total of 25,000 individuals, following the model $y_{ref} = Ma_{ref} + e_{ref} = g_{ref} + e_{ref}$, where M is a matrix of which element $M_{ij}$ corresponds to the centered genotype of individual i at QTL j, $a_{ref} = [a_{ref,1} ... a_{ref,q}]$ is the vector of the q = 512 additive QTL effects. The QTL effects $a_{ref,j}$ were initially sampled for every j = 1, ..., q, as independent and identically distributed (i.i.d.) from a normal distribution $\alpha_{ref,j} \sim N(0, \sigma^2)$, with

$\sigma^2 = \frac{\sigma_{y_{ref}}^2 \times h^2}{\sum_{j=1}^q 2p_j(1-p_j)}$, such that h² is the heritability,

$\sigma^2$ is the phenotypic variance, and p is the allele frequency. The vector of the true genetic values (TGV) was computed as $g_{ref} = Ma_{ref}$, with $g_{ref} \sim N(0, \sigma^2)$ and $\sigma^2$ as the genetic variance; $e_{ref} \sim N(0, \sigma^2)$ is the vector of random errors where $\sigma^2$ was chosen to ensure that the $\sigma^2 = \sigma^2 + \sigma^2$. The genotypes and the reference trait were simulated using the GenEval package (version 1.0.0, https://github.com/bcuyabano/GenEval; Cuyabano, 2021) implemented in R (https://www.r-project.org/).

Five dependent traits were simulated conditional to the TGV of the reference trait as $y_t = g_t + e_t$, for every t = 1, ..., 5, such that $g_t = f_t(g_{ref}) + g_{t-2}$, where $f_t$ is the function describing the (potentially) nonlinear relationship between the TGV of 2 traits; $g_{t-2}$ is the vector of the GV specific to each dependent trait, simulated as $g_{t-2} = M_t a_t$, where $M_t$ is the genotype matrix of QTLs specific to trait t and different from the common QTLs and different between the dependent traits; and $\alpha_t \sim N(0, \sigma^2)$ is the vector of the corresponding additive QTL effects. $g_t \sim N(0, \sigma^2)$ was normalized as $\sigma^2 = \sigma^2 \times h^2$, where h² is the heritability. The error vector was simulated as $e_t \sim N(0, \sigma^2)$.

For the simulation of linearly correlated traits, the genetic and error (co)variance matrices were $[\sigma^2_{g_{ref}} \sigma^2_{g_{ref,2}} \sigma^2_{g_{ref,3}} \sigma^2_{g_{ref,4}} \sigma^2_{g_{ref,5}}] \otimes G$ and $[\sigma^2_{e_{ref}} \sigma^2_{e_{ref,1}} \sigma^2_{e_{ref,2}} \sigma^2_{e_{ref,3}} \sigma^2_{e_{ref,4}} \sigma^2_{e_{ref,5}}] \otimes I$, respectively, where $\otimes$ is the Kronecker product. We fixed $\sigma^2 = 20$, without loss of generality for the proof-of-concept results on simulated data, to allow a fair comparison between traits. Different values of h² (0.6, 0.3, 0.15, 0.05) and different numbers of specific QTLs (0, 10, 50, 250) were considered for the traits. Although different levels of h² and number of specific QTL were tested, the values were the same between traits of each set. Four relationships were considered for the $f_t$:

1. Linear: $f_t(g_{ref}) = 1$

$1$

$1$

4. Sinusoidal: $f_t(g_{ref}) = \sin(p_1 g_{ref})$

The coefficients used were $p_1 = 1$, $p_2 = 1$, and $p_3 = 2$. To assess the repeatability of our findings, we simulated 50 replicas of the complete dataset under each combination of different levels of h² and number of QTL.

To train and evaluate the DL models, the dataset was split into 3 sets: training (80%), validation (10%), and test sets (10%). The test set consisted of individuals for whom PGV was to be obtained in the absence of phenotypic records, whereas both the training and validation sets contained individuals with genotypes and phenotypic records. Different from statistical models, DL requires an internal validation set with the complete information to fine-tune the model parameters. For the GBLUP model, the training and validation sets were combined to fit the
model, and the same test set (10% of individuals) was kept for evaluation.

## Real Data

The real dataset consisted of 113,599 genotyped Holstein cows, with a total of 53,469 SNPs after editing (original SNP data from the EuroGMD v1, a customized Illumina genotyping microarray manufactured by Illumina Inc. that contains ~70,000 SNPs). For the performance records, we used the yield deviations (YD) of 33 traits, all of them available for every one of the 113,599 cows. These included 5 dairy production traits (heritability 38%–63%), SCS (heritability 19%), 6 fertility traits (heritability 1%–7%), and 21 dairy morphology traits (heritability 10%–51%). The YD consisted of the phenotype corrected for every nongenetic effect from single-trait models used in the French genetic evaluation (effects such as the contemporary group, lactation number, and age of the cow, among others, were adequately considered according to each one of the traits; for details of the French genetic evaluation, see GenEval, 2024) and were obtained using the single-step GBLUP (ssGBLUP). For the real dataset, we split the data into 3 subsets (training [82%], validation [9%], and test [9%]) based on the generation structure of the animals. To achieve this, we grouped animals according to their birth years, assigning older generations to the training set, intermediate generations to the validation set, and the most recent generations to the test set. This chronological separation ensures that the model is evaluated on progressively newer animals, reflecting real-world genetic evaluation scenarios. We ensured that no cow had paternal half-sisters in different subsets to prevent information leakage.

## Multitrait Genomic Prediction

We used 3 models to perform GP: a baseline GBLUP, a DL model, and a hybrid DLGBLUP model that combines the 2 previous ones.

**Genomic Best Linear Unbiased Prediction.** The GB-LUP model is one of the most popular statistical methods used to predict the GV of genotyped individuals using their genomic relationship matrix. The model considered for multitrait GBLUP is as follows:

$$ \mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{g} + \mathbf{e}, \quad [1] $$

where $\mathbf{Y} = [\mathbf{Y}_1 \dots \mathbf{Y}_T]$ is a matrix of the phenotype vectors for all trait included in the model, of dimension $n \times T$, where $n$ is the number of records and $T$ the number of traits; $\mathbf{X}$ is a design matrix for the fixed effects, of dimension $n \times p$, where $p$ is the number of fixed effects,

$\boldsymbol{\beta} = [\boldsymbol{\beta}_1 \dots \boldsymbol{\beta}_t]$, such that $\boldsymbol{\beta}_t$ is the vector of fixed effects for trait $t = 1, \dots, T$; $\mathbf{g} = [\mathbf{g}'_1 \dots \mathbf{g}'_T]'$ is a matrix of genetic values of dimension $n \times T$, such that $\mathbf{g}_t$ is the vector of genetic values for trait $t$, with distribution $\mathbf{g} \sim N(0, \mathbf{\Sigma}_g \otimes \mathbf{G})$, in which $\mathbf{\Sigma}_g$ is the additive genetic (co) variance matrix of the traits; $\mathbf{G}$ is the genomic relationship matrix (VanRaden, 2008); $\otimes$ represents their Kronecker product; and $\mathbf{e} = [\mathbf{e}'_1 \dots \mathbf{e}'_T]'$ is the matrix of normally distributed random errors with $\mathbf{e} \sim N(0, \mathbf{R} \otimes \mathbf{I})$, in which $\mathbf{R}$ is the error (co)variance matrix between the traits. We used the R (4.1.1) package BGLR (1.0.8; Pérez-Rodríguez and de los Campos, 2022) to fit the multitrait GBLUP on the simulated data, but here we did not consider a fixed effect, so $\boldsymbol{\beta}$ is only the vector of the intercepts. We ran the models with 25,000 iterations, discarding the first 5,000 samples as burn-in with a thin 10. For the real dataset, a single-trait ssGBLUP analysis was performed using the HSSGBLUP software (Tribout et al., 2020), developed for the French cattle genetic evaluation.

**Deep Learning Model.** The multilayer perceptron (MLP; Goodfellow et al., 2016) is a feedforward network that maps an input to an output using learned parameters. The model transforms the input by passing through multiple connected layers (hidden layers) up to the output. In a general execution, each node is characterized by the following equation:

$$ y_j = f\left(\sum_{i=1}^{I} w_{j,i}x_i + b_j\right), \quad [2] $$

where $y_j$ is the scalar output of node j, $b_j$ is the bias, $x_i$ represents the output of node i in $\{1, \dots, I\}$ of the previous layer, $w_{j,i}$ is the weight assigned to this $x_i$, and f is a nonlinear activation function.

The proposed DL model, illustrated in Figure 1, comprises 2 distinct modules: SNPs2Trait and Trait2Trait. The first module, SNPs2Trait, takes the SNP data as input and provides a first predicted GV ($\text{PGV}_1$) for all traits as an output by capturing the relationship between the genomic data and the phenotype, defined as

$$ \text{PGV}_1 = f_{\text{MLP}}(\mathbf{X}), $$

where $\mathbf{X}$ represents the input SNP matrix. The hidden layers enable the construction of a latent genome representation common to all traits. This latent representation is used in the output layer to predict each trait separately. The second module, Trait2Trait, takes as input the $\text{PGV}_1$ of all traits and outputs a new predicted GV ($\text{PGV}_2$) for all traits by capturing the relationship between all traits that will be considered as the final prediction of
<figure><img src="image_3.png" /><figcaption>Figure 1. Deep learning genomic prediction frameworks. Pheno = phenotype.</figcaption></figure>

the model. The Trait2Trait module refines PGV<sub>-1</sub> through additional transformations:

$$PGV_{-2} = g_{\text{MLP}}(PGV_{-1}),$$

where $g_{\text{MLP}}$ represents the neural network in this module.

During training, the optimization process adjusts the weights and biases of the MLP by minimizing a loss function that quantifies the difference between the predicted and true outputs on the training set. The mean squared error (MSE) and mean absolute error (MAE) losses are commonly used for regression, whereas cross entropy loss is used for classification. Here, we used Huber loss (Hastie et al., 2001) due to its combination of the advantages of both MSE and MAE, providing a balance between sensitivity to small errors and robustness to outliers. This loss is described as

$$L(z) = \begin{cases} \frac{1}{2}(z)^2 & \text{for } |z| \le \delta, \\ |z| - \frac{1}{2}, & \text{otherwise.} \end{cases}, \quad [3]$$

where z is the error and δ (fixed to 1) is a threshold parameter that dictates the transition between the 2 losses. In practice, stochastic gradient descent is used to update the parameters (weights and biases). We used Adam optimizer (adaptive moment estimation) with a learning rate of 10<sup>-4</sup>. Two training approaches were used. The first approach involved training the 2 modules as one model so the loss function was the sum of loss from each module:

$$Loss_{-2} = \sum_{t} L(y_{t} - PGV_{-1}) + L(y_{t} - PGV_{-2}), \quad [4]$$

where $y_t$ is the true phenotype value for trait t. Incorporating both $L(y_t - PGV_{-1})$ and $L(y_t - PGV_{-2})$ ensures that the model learns meaningful intermediate representations—specifically, the genetic value—rather than merely transforming the genotype into a phenotype prediction through a larger network. If only $L(y_t - PGV_{-2})$ were used, the model would resemble the SNP2Trait module, treating the input genotype as a direct predictor of the phenotype without explicitly capturing the GV representation. The second approach involved training modules sequentially in terms of gradient flow: the first module was trained initially using the loss function $L(y_t - PGV_{-1})$, and once this was completed, its parameters were fixed by stopping the gradient flow. This procedure isolated it from the training phase of the second module using $L(y_t - PGV_{-2})$, maintaining the integrity of the initial module's parameters as the later module was trained.

The SNPs2Trait module consists of one hidden layer (with 400 neurons), and the Trait2Trait module consists of 2 hidden linear layers (with 400 and 256 neurons, respectively). The LeakyReLU activation function was used with a fixed negative slope of 0.1. The model was trained for 100 maximum epochs, with each epoch consisting of a single pass through the training dataset, using a batch size of 200. An early stopping was employed to terminate training if there was no improvement in validation loss after 10 epochs. The choice of optimizer, learning rate, activation function and number of neurons in hidden layers was made after a grid search (Supplemental Table S1, see Notes). The model was implemented using pytorch (1.10.2) in Python and trained on a single GPU with 48 GB memory (NVIDIA A40).

**DLGBLUP Model.** We proposed a hybrid version of the previous model called DLGBLUP (Figure 2), where the first module is replaced with GBLUP. The reason behind this variation is that the first module, SNPs2Trait, performs a similar function to GBLUP by learning the patterns between genotype and phenotype to predict the GV. We wanted to test the effect of starting with the results of GBLUP to make the final predictions using the Trait2Trait module, which remains unchanged.

**Intertrait Relationship Representation.** In order to represent graphically the relationship between the reference trait and a dependent trait, a MeanTrait2Trait module was proposed (Figure 3). This module outputs predictions.

<figure><img src="image_4.png" /><figcaption>Figure 2. The DLGBLUP genomic prediction frameworks. The GBLUP framework is shown in red, and the DL framework is shown in green.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3. Schematic representation of the MeanTrait2Trait model. The model takes the predicted genetic values for the reference trait (PGVref) as input to predict the PGV of the dependent traits (PGVT_3). The plot on the right illustrates the relationship between PGVref and a dependent trait (PGV2_3), demonstrating a nonlinear trend.</figcaption></figure>

tions of the dependent traits which we call PGV₃ by taking PGV of reference trait (PGVref) as input. As such, MeanTrait2Trait can capture the mean relationship between reference trait (input) and the other dependent traits. The training of MeanTrait2Trait was done using $Loss_3 = \sum_t L(PGV_t - PGV_{t-3})$. This last module is used for illustration purposes only.

## Impact of QTL Exclusion and Noncausative SNP Addition

We first considered the ideal case in which only the true QTL were provided to the prediction methods, serving as our baseline. To evaluate the impact of including SNPs that were not QTL on prediction accuracy, we gradually added noncausative SNPs to the input genomic data. This resulted in scenarios where the genomic composition included 10%, 25%, 50%, and 75% non-QTL SNPs. Lastly, a final scenario consisted of using only noncausative SNPs, excluding all the QTL. A comparison of the results from a model excluding noncausative mutations, but including only half of the QTL, was also included in this study. These scenarios offer a clear insight into how the signal-to-noise ratio affects the prediction of PGV and the effect of relationships between traits.

## Genomic Selection

After performing the multitrait GP, we computed a selection index (SI) as follows:

$$SI = \sum_{t=1}^{T} \left( w_t + \sum_{j=1, j \neq t}^{T} \rho_{t,j} w_j \right) g_t, \quad [5]$$

in which $w_t$ are the weights given to traits $t = 1, \dots, T$; $\rho_{t,j}$ is the empirical correlation between the GV of traits $t$ and $j$; and $g_t$ are the GV vector of trait $t = 1, \dots, T$. For this study, we assigned equal weights to all traits in order to compare the general evolution of the genetic gain

based on the additive PGV obtained with the different models. Based on the SI, at each generation the top 10% eligible males were selected and mated with all females from the latest generations, in order to generate the next generation, maintaining at each new generation a female-to-male proportion of 80% to 20%. Males were kept for breeding across 4 generations, whereas females were kept for 2 generations only. This selection and mating design was chosen to mimic the breeding program of a dairy cattle population. A base population of independent individuals (i.e. generation zero, or $G_0$) was considered to start the selection process, followed by 7 simulated generations under selection referred to as $G_1$ to $G_7$.

This simulation scheme under selection allowed us to both study the evolution of the nonlinear genetic relationships between traits when the population is under selection, and compare the genetic gain of selection, based on the additive GV obtained with the different models, using their values to construct the SI in Equation 5. Because of the way we simulate the data, the PGV obtained by all prediction models could express nonadditive effects. Because we are interested in the additive PGV, we computed it as follows: at each generation we simulated the genotype of 20 offspring for each eligible male based on random mating schemes with the females. The PGV of these simulated offspring were then predicted with the different models (GBLUP, DL, and DLGBLUP) trained on a dataset where all individuals until the latest generation were genotyped and phenotyped. The additive PGV of the eligible males were finally computed as the average of the PGV of their simulated offspring.

The TGV for the reference traits on all subsequent generations were simulated using the same original QTL effects of $G_0$, and the other dependent traits were simulated as previously described in the “Simulated Data” section. The number of individuals in each of the generations $G_1$ to $G_7$ was maintained at 5,000. With each new generation, the reference population to perform the prediction models was updated to comprise all individuals from $G_0$ to the latest with both genotypes and phenotypic records, and then used to perform the genetic evaluation. For DL and DLGBLUP, the genetic evaluation model with a new generation was initialized using the weights (DL parameters) obtained at the evaluation of the previous generation.

## Evaluation Metrics of the Predicted Genetic Values of Individual Traits

To evaluate the performance of each model in predicting GV, we used the Pearson correlation coefficient to calculate the prediction accuracy for each individual trait, as well as the MSE between the TGV and PGV, and performed a visual assessment of the relationship
form between traits. We also computed how much, in
terms of variance percentage, a deviation from linear-
ity represented, compared with the variation in the PGV
of dependent traits. This was done by calculating the
variance of the linear or nonlinear predictions using the
MeanTrait2Trait model of the dependent trait and divid-
ing it by the variance of the predicted genetic value of
the same trait multiplied by 100. Although this variance
percentual is not strictly a formal manner to decompose
the variance as would be a sequential ANOVA, it does
give an approximation of how much linear variance can
be explained by the nonlinear relationship between traits.
Although a regression of the nonlinear function would be
the ideal solution, the MeanTrait2Trait module does not
output a formal function to describe the nonlinear rela-
tionship between the 2 traits, but only a mean per range
of values, allowing a visual assessment of nonlinearity.
Therefore, we opted to calculate the variance percentual,
which although not a precise measure, provides an ap-
proximation of the proportion of variance explained. To
evaluate the performance of each model with respect to
genomic selection, we used the genetic gain computed
as the difference between the mean of the TGV for indi-
viduals in G₀ and the mean of the TGV for individuals in
the next generation.

RESULTS

Here we present the results of DL-based model with
MLP architecture and trained following the second train-
ing approach (sequential training) which gave better
results. The reported results are based on the testing set.
For simplicity, we named the traits based on the form of
relationship they have with the reference trait; for ex-
ample, quadratic traits correspond to the trait with f as
quadratic. We used the phenotype set with heritability 0.3
and 50 specific QTLs as default. The results of sinusoidal
trait can be found in Supplemental Figure S1 (see Notes),
as it was used to test an extreme case of nonlinearity.

Genomic Prediction

*Prediction Accuracy.* Figure 4 presents the perfor-
mance comparison of the GBLUP, DL, and DLGBLUP
models using only the QTL as input data, with respect to
the prediction accuracy and the MSE, over the 50 repli-
cates performed.

At the additive effect level, GBLUP demonstrated
superior prediction accuracy compared with the DL
SNP2Trait model for reference and linear traits. The
opposite was shown for the quadratic traits. For logistic

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Performance comparison of different prediction models for various traits. GBLUP (red), deep learning (DL; yellow); with the 2 modules SNP2Trait and Trait2Trait and DLGBLUP (green). Error bars represent SD across replicates. The prediction was made using genomic data with only QTL.</figcaption></figure>
traits, all methods exhibited similar results. Nonetheless, DL SNPs2Trait presented a consistently lower MSE for all traits, except the reference trait, as expected due to the shrinkage of the marker effects in the GBLUP procedure.

Independent of the procedure used to obtain PGV_1 (either GBLUP or DL SNP2Trait), the application of the Trait2Trait module always improved the prediction accuracy of PGV_2 (Figure 4). The prediction accuracy increased by 0.01 to 0.07 for DL Trait2Trait module, and by 0.01 to 0.2, depending on the type of nonlinear relationships between traits. After passing through the Trait-2Trait module, no further improvement was provided for the prediction accuracy of the linear traits (reference included), which were already well predicted specially by GBLUP.

We found that DLGBLUP outperformed DL Trait2Trait by a margin of 0.01 to 0.02 in terms of prediction accuracy. Of all the evaluated models, DLGBLUP Trait2Trait consistently yielded the lowest MSE across all dependent traits, with a more reliable and precise estimation and a prediction accuracy equal or higher than all the other models.

**Intertrait Relationship Representation.** The Mean-Trait2Trait module captured the relationship between traits correctly (Figure 5), whether it was linear or nonlinear, using the PGV<sub>ref</sub> predicted from any prediction method. Moreover, the shape of the PGV was adjusted according to the trait's relationship through the Trait-2Trait part. The SNPs2Trait module was able to identify a slight nonlinearity between the traits but not the complete shape detected by Trait2Trait. A more accurate input (i.e., the GBLUP predictions) to the Trait2Trait module enabled the prediction of a clearer and more precise relationship between traits, as illustrated in Figure 5 and Supplemental Figure S2 (see Notes). We observed

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison plots of predicted genetic values for quadratic trait with SNP2Trait and Trait2Trait modules of the DL model and GBLUP and Trait2Trait modules of DLGBLUP. The true GV is shown in blue, DL SNP2Trait prediction in light yellow, DL Trait2Trait prediction in orange, GBLUP prediction in red, and Trait2Trait DLGBLUP in green.</figcaption></figure>

<table><caption>**Table 1.** Variance percentual of linear and nonlinear relationship prediction for simulated data</caption><thead><tr><th rowspan="2">Trait</th><th colspan="2">Variance %</th></tr><tr><th>Linear</th><th>Nonlinear</th></tr></thead><tbody><tr><td>Linear</td><td>24.4</td><td>21.58</td></tr><tr><td>Quadratic</td><td>0.26</td><td>34.68</td></tr><tr><td>Logistic</td><td>17.18</td><td>24.17</td></tr></tbody></table>

that GBLUP transformed the relationship to be linear, sometimes identifying a level of linear relationship (logistic trait) or completely missing any relationship between traits (quadratic trait). In contrast, Trait2Trait was able to detect the genetic relationships between traits. When examining the contribution of this relationship to the total GV, we noticed a high percentage of nonlinear variance contribution for nonlinear traits, with a large difference between the nonlinear and linear contributions to this variance, especially for quadratic traits. There is an important contribution from the linear relationship for linear models predicting linear traits, which is higher than the nonlinear contribution (Table 1).

**Effect of Heritability, Number of Specific QTL, and Level of Genetic Relationship.** Heritability directly influences the quality of GV predictions from additive effects. The higher the heritability, the higher the prediction accuracy of GBLUP, as well as the prediction accuracy made by the Trait2Trait part and its improvement in comparison to GBLUP for the nonlinear traits, reaching 0.17 for the quadratic traits with a heritability of 0.6 (Table 2).

Moreover, we observed an effect of the number of specific QTL that control the genetic relationship (Figure 6). The higher the number of specific QTL, the less clear the relationship between traits. A higher genetic relationship (a lower number of specific QTL in our case) intensifies the nonlinear effects, reducing the prediction accuracy of GBLUP while simultaneously increasing the potential for improvement via the Trait2Trait component, as shown in Table 3.

## Impact of QTL Exclusion and Noncausative SNP Addition

Table 4 presents the impact on the prediction accuracy of incorporating SNPs that were not QTL in the genomic data used. The most accurate predictions were achieved when the input data comprised exclusively all the QTLs implicated in the TGV of a trait, as expected, with DLGBLUP achieving the best predictions. When excluding half of the QTLs or when adding noncausative SNPs to the input genotypes, the prediction accuracy decreased. The decrease in prediction accuracy when noncausative SNPs
<table><caption>**Table 2.** Prediction accuracy comparison of different prediction models with different heritability with 50 specific QTL</caption><thead><tr><th rowspan="2">Heritability</th><th rowspan="2">Model</th><th colspan="4">Trait</th></tr><tr><th>Referent</th><th>Linear</th><th>Quadratic</th><th>Logistic</th></tr></thead><tbody><tr><td rowspan="2">0.05</td><td>GBLUP</td><td>0.83</td><td>0.8</td><td>0.8</td><td>0.84</td></tr><tr><td>DL</td><td>0.72</td><td>0.67</td><td>0.7</td><td>0.74</td></tr><tr><td rowspan="2">0.15</td><td>DLGBLUP</td><td>0.82</td><td>0.79</td><td>0.8</td><td>0.83</td></tr><tr><td>GBLUP</td><td>0.93</td><td>0.94</td><td>0.82</td><td>0.93</td></tr><tr><td rowspan="2">0.3</td><td>DL</td><td>0.88</td><td>0.89</td><td>0.81</td><td>0.89</td></tr><tr><td>DLGBLUP</td><td>0.93</td><td>0.93</td><td>0.87</td><td>0.92</td></tr><tr><td rowspan="2">0.6</td><td>GBLUP</td><td>0.97</td><td>0.97</td><td>0.84</td><td>0.96</td></tr><tr><td>DL</td><td>0.95</td><td>0.95</td><td>0.91</td><td>0.94</td></tr><tr><td rowspan="2"></td><td>DLGBLUP</td><td>0.97</td><td>0.97</td><td>0.94</td><td>0.95</td></tr><tr><td>GBLUP</td><td>0.99</td><td>0.99</td><td>0.79</td><td>0.96</td></tr><tr><td rowspan="2"></td><td>DL</td><td>0.98</td><td>0.98</td><td>0.96</td><td>0.97</td></tr><tr><td>DLGBLUP</td><td>0.99</td><td>0.99</td><td>0.96</td><td>0.98</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison plots of the true genetic relationship between traits with different number of specific QTLs. The reference trait is represented in the x-axis, and all dependent traits in the y-axis. The blue points represent the true GV of a population.</figcaption></figure>
<table><caption>**Table 3.** Prediction accuracy comparison of different prediction models with different number of specific QTL and heritability 0.3</caption><thead><tr><th rowspan="2">QTLs (n)</th><th rowspan="2">Model</th><th colspan="4">Trait</th></tr><tr><th>Referent</th><th>Linear</th><th>Quadratic</th><th>Logistic</th></tr></thead><tbody><tr><td rowspan="3">0</td><td>GBLUP</td><td>0.99</td><td>0.99</td><td>0.13</td><td>0.94</td></tr><tr><td>DL</td><td>0.97</td><td>0.97</td><td>0.85</td><td>0.93</td></tr><tr><td>DLGBLUP</td><td>0.99</td><td>0.99</td><td>0.95</td><td>0.97</td></tr><tr><td rowspan="3">10</td><td>GBLUP</td><td>0.98</td><td>0.98</td><td>0.42</td><td>0.95</td></tr><tr><td>DL</td><td>0.96</td><td>0.96</td><td>0.82</td><td>0.94</td></tr><tr><td>DLGBLUP</td><td>0.98</td><td>0.97</td><td>0.9</td><td>0.96</td></tr><tr><td rowspan="3">50</td><td>GBLUP</td><td>0.97</td><td>0.97</td><td>0.84</td><td>0.96</td></tr><tr><td>DL</td><td>0.95</td><td>0.95</td><td>0.91</td><td>0.94</td></tr><tr><td>DLGBLUP</td><td>0.97</td><td>0.97</td><td>0.94</td><td>0.95</td></tr><tr><td rowspan="3">250</td><td>GBLUP</td><td>0.92</td><td>0.93</td><td>0.9</td><td>0.93</td></tr><tr><td>DL</td><td>0.89</td><td>0.89</td><td>0.86</td><td>0.89</td></tr><tr><td>DLGBLUP</td><td>0.92</td><td>0.92</td><td>0.91</td><td>0.93</td></tr></tbody></table>

were introduced was progressive as more and more SNPs were included to the genomic data, resulting in a decrease in the advantage of DLGBLUP. Excluding all QTL and keeping only noncausative SNPs in LD with the QTL resulted in PGV with accuracies ranging from 0.55 to 0.7 and an improvement made by DLGBLUP ranging from 0.01 to 0.03, when compared with both DL and GBLUP.

## Genomic Selection

**Genetic Relationship Across Generations.** The prediction accuracy of the predicted additive genetic effects for all methods was similar, ranging between 0.9 and 0.97,

as shown in Table 5. Figure 7A shows the evolution of the genetic relationships between the TGV of the traits over the 7 generations, under a multitrait selection using the SI based on the additive GV obtained with pure DL, using only the QTLs as input data. As expected, the relationship between the reference trait and the linear trait remained linear over all the generations under selection. In contrast, the relationship between the reference trait and the quadratic trait evolved from nonlinear to linear very quickly, with its nonlinearity being almost imperceptible from $G_1$. The logistic trait remained nonlinear, although the shape of its relationship with the reference trait changed completely with respect to the original one

<table><caption>**Table 4.** Prediction accuracy comparison of different prediction models using genomic data with different numbers of SNPs, a $h^2$ of 0.3 and 50 specific QTLs</caption><thead><tr><th rowspan="2">Item</th><th rowspan="2">Method</th><th colspan="4">Trait</th></tr><tr><th>Referent</th><th>Linear</th><th>Quadratic</th><th>Logistic</th></tr></thead><tbody><tr><td rowspan="3">Half QTL</td><td>GBLUP</td><td>0.82</td><td>0.76</td><td>0.66</td><td>0.78</td></tr><tr><td>DL</td><td>0.82</td><td>0.75</td><td>0.72</td><td>0.77</td></tr><tr><td>DLGBLUP</td><td>0.81</td><td>0.74</td><td>0.72</td><td>0.77</td></tr><tr><td rowspan="3">Just QTL</td><td>GBLUP</td><td>0.97</td><td>0.97</td><td>0.76</td><td>0.96</td></tr><tr><td>DL</td><td>0.95</td><td>0.95</td><td>0.95</td><td>0.97</td></tr><tr><td>DLGBLUP</td><td>0.97</td><td>0.97</td><td>0.96</td><td>0.98</td></tr><tr><td rowspan="3">10% SNPs</td><td>GBLUP</td><td>0.95</td><td>0.95</td><td>0.8</td><td>0.95</td></tr><tr><td>DL</td><td>0.94</td><td>0.94</td><td>0.9</td><td>0.94</td></tr><tr><td>DLGBLUP</td><td>0.94</td><td>0.94</td><td>0.91</td><td>0.94</td></tr><tr><td rowspan="3">25% SNPs</td><td>GBLUP</td><td>0.91</td><td>0.87</td><td>0.75</td><td>0.9</td></tr><tr><td>DL</td><td>0.87</td><td>0.89</td><td>0.82</td><td>0.88</td></tr><tr><td>DLGBLUP</td><td>0.9</td><td>0.9</td><td>0.84</td><td>0.89</td></tr><tr><td rowspan="3">50% SNPs</td><td>GBLUP</td><td>0.86</td><td>0.86</td><td>0.7</td><td>0.86</td></tr><tr><td>DL</td><td>0.82</td><td>0.83</td><td>0.74</td><td>0.82</td></tr><tr><td>DLGBLUP</td><td>0.85</td><td>0.85</td><td>0.78</td><td>0.85</td></tr><tr><td rowspan="3">75% SNPs</td><td>GBLUP</td><td>0.84</td><td>0.84</td><td>0.68</td><td>0.84</td></tr><tr><td>DL</td><td>0.78</td><td>0.79</td><td>0.71</td><td>0.79</td></tr><tr><td>DLGBLUP</td><td>0.82</td><td>0.83</td><td>0.75</td><td>0.82</td></tr><tr><td rowspan="3">All SNPs</td><td>GBLUP</td><td>0.82</td><td>0.82</td><td>0.67</td><td>0.82</td></tr><tr><td>DL</td><td>0.76</td><td>0.77</td><td>0.69</td><td>0.78</td></tr><tr><td>DLGBLUP</td><td>0.8</td><td>0.81</td><td>0.73</td><td>0.81</td></tr><tr><td rowspan="3">No QTL</td><td>GBLUP</td><td>0.63</td><td>0.67</td><td>0.56</td><td>0.7</td></tr><tr><td>DL</td><td>0.61</td><td>0.64</td><td>0.57</td><td>0.67</td></tr><tr><td>DLGBLUP</td><td>0.63</td><td>0.67</td><td>0.59</td><td>0.69</td></tr></tbody></table>
<table><caption>**Table 5.** Prediction accuracy of additive genetic value comparison of different prediction models</caption><thead><tr><th>Item</th><th>GBLUP</th><th>DL</th><th>DLGBLUP</th></tr></thead><tbody><tr><td>Reference</td><td>0.97</td><td>0.95</td><td>0.97</td></tr><tr><td>Linear</td><td>0.97</td><td>0.96</td><td>0.97</td></tr><tr><td>Quadratic</td><td>0.9</td><td>0.92</td><td>0.93</td></tr><tr><td>Logistic</td><td>0.96</td><td>0.95</td><td>0.96</td></tr></tbody></table>

in $G_0$. These evolutions of the shape of the relationships persisted whether the selection was based on any of the predicted additive GV (DL, GBLUP, or DLGBLUP; Supplemental Figure S3, see Notes).

**Genetic Gain.** Genetic gain was achieved for all the traits, following the selection of the top 10% males (Figure 7B) based on the SI with any of the additive PGV. The magnitude of this gain varied depending on the form of relationship, with the most substantial gains observed for the logistic trait. For both reference and

linear traits, the progress based on the SI using the additive GV from all models were equivalent. Conversely, for each of the nonlinear traits, the genetic progress achieved based on the SI using the additive GV from either DL or DLGBLUP, or both for some traits, was greater than the genetic progress achieved using the additive GV from GBLUP.

## Application on Real Data

The DLGBLUP method was applied on real data. The YD was used instead of phenotype in simulated data, and the results of unitrait GBLUP of 33 traits were used as input for Trait2Trait module. The results of 6 traits are presented: milk yield, fat yield, and protein content as production traits; SCS as a health trait; conception rate of cows as a fertility trait; and depth of furrow as a morphology trait (for other traits, refer to Supplemental Table S1). On test set, Trait2Trait yielded a lower MSE

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparative analysis of traits progression across generations. (A) Plots of the relationship between the TGV of reference traits (x-axis) and dependent traits (y-axis) over 8 generations based on DL predictions. (B) Comparison between the genetic gain of 7 generations with 15% male selection using the true and the predicted additive GV with GBLUP, DL, and DLGBLUP for all traits.</figcaption></figure>
<table><caption>**Table 6.** Prediction accuracy and MSE comparison between GBLUP and Trait2Trait module of DLGBLUP model for multiple traits for real data on test set</caption><thead><tr><th rowspan="2">Trait</th><th rowspan="2">Heritability</th><th colspan="2">Prediction accuracy</th><th colspan="2">MSE</th></tr><tr><th>GBLUP</th><th>Trait2Trait</th><th>GBLUP</th><th>Trait2Trait</th></tr></thead><tbody><tr><td>Milk yield</td><td>0.38</td><td>0.65</td><td>0.65</td><td>0.65</td><td>0.59</td></tr><tr><td>Fat yield</td><td>0.38</td><td>0.62</td><td>0.62</td><td>0.8</td><td>0.64</td></tr><tr><td>Protein content</td><td>0.63</td><td>0.89</td><td>0.89</td><td>0.21</td><td>0.21</td></tr><tr><td>SCC</td><td>0.19</td><td>0.34</td><td>0.35</td><td>1.38</td><td>0.92</td></tr><tr><td>Conception rate cows</td><td>0.04</td><td>0.12</td><td>0.12</td><td>1.85</td><td>1.05</td></tr><tr><td>Depth of furrow</td><td>0.23</td><td>0.38</td><td>0.38</td><td>1.23</td><td>0.88</td></tr></tbody></table>

for all traits; however, no improvement was observed in terms of prediction accuracy (Table 6 and Supplemental Table S2, see Notes).

Even if the relationship seemed to be linear when plotted over the PGV on population level, MeanTrait2Trait was able of detecting a form of nonlinearity between different traits (Figure 8). For example, the relationship between milk yield and fat yield demonstrates a positive, linear-like trend, whereas the relationships between the other traits exhibit a more nonlinear pattern. However, when considering the contribution percentage, both linear and nonlinear relationships contribute only a small percentage (less than 1%) for most traits, except for milk yield (~6%). The contribution of nonlinearity was higher than that of linearity for all traits except for the conception rate, where the contribution of linearity was twice as large as the contribution of nonlinearity (Table 7).

## DISCUSSION

In this study, we proposed the use of DL to predict potential nonlinear genetic relationships between traits. To do so, we developed DLGBLUP, a novel hybrid model that integrates deep learning and statistical methods to enhance multitrait prediction by capturing potential nonlinear genetic relationships between traits. The DLGBLUP model was able to predict nonlinear relationships in both simulated and real data from the French Holstein dairy cattle population. With simulated data, we have shown that selection based on the additive PGV from DLGBLUP achieved a greater genetic progress than selection based on additive PGV from GBLUP.

Nonlinear genetic relationships between traits have been largely overlooked in GP, both in DL and statistical methods. Deep learning has mainly been explored as an alternative to statistical models, yet its performance on real datasets has not substantially surpassed GBLUP (Montesinos-López et al., 2021). One reason could be that GP aims to predict the additive GV—unobserved in phenotype data—without prior assumptions for DL models. Unlike other fields, GP requires DL models to infer genotype-phenotype relationships through hidden layers, often lacking biological interpretability. Some studies have addressed this by integrating biological insights or combining DL with statistical methods, leading to improved performance (van Hilten et al., 2021; Lee et al., 2023).

By accounting for the possibility of nonlinear relationships between traits, we have shown with our study that the DL and DLGBLUP models outperformed GBLUP on the simulated nonlinear traits, successfully capturing the shape of the relationships without any prior assumptions about their form, especially if the provided input is more accurate. In addition, we illustrated the critical role of an accurate SNPs selection step in minimizing input noise without eliminating important information. However, the model's performance depends on the distribution of SNPs and QTLs in the genomic data. This limitation arises not from the model's ability to identify relationships, but from the quality of input data. If PGV<sub>i</sub> is noisy or poorly predicted due to genotype quality, the Trait2Trait module will have difficulty accurately capturing genetic relationships, leading to decreased prediction accuracy and reliability. This underscores the importance of high-quality genomic information for optimizing model performance. A notable advantage of DLGBLUP is its complementarity with other single-trait or multitrait statistical methods, such as GBLUP or even advanced approaches like the single-step method, which can predict values for nongenotyped individuals using pedigree information. When relationships between traits are highly nonlinear, this complementarity enables DLGBLUP to achieve superior results, including higher prediction accuracy, lower prediction error, and better deciphering of

<table><caption>**Table 7.** Contribution of linear and nonlinear relationship prediction for real data</caption><thead><tr><th rowspan="2">Trait</th><th colspan="2">Contribution %</th></tr><tr><th>Linear</th><th>Nonlinear</th></tr></thead><tbody><tr><td>Milk yield</td><td>5.68</td><td>6.39</td></tr><tr><td>SCC</td><td>0.16</td><td>0.27</td></tr><tr><td>Conception rate cows</td><td>0.33</td><td>0.15</td></tr><tr><td>Depth of furrow</td><td>0.13</td><td>0.59</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Plots of the predicted genetic values with GBLUP, DLGBLUP and the relationship representation between traits on real data. Red dots represent the predicted genetic value using GBLUP and green dots represent the predicted genetic value using DLGBLUP. Black lines indicate the predicted relationship representation using MeanTrait2Trait, and blue lines indicate the linear relationship between 2 traits predicted using GBLUP. The traits in (y-axis) represents the dependent trait predicted from the trait in x-axis (input) for MeanTrait2Trait model.</figcaption></figure>

the relationships between traits, compared with using DL alone. The DLGBLUP method can process the outputs of routine genetic evaluations by postprocessing the breeding values of multiple traits. It leverages the phenotypes, genotypes, and initial breeding values to identify nonlinear genetic relationships through its Trait2Trait module, which detects and accounts for these complex interactions between traits.

Another important feature of the new approach is that it comes with no substantial computational overload. Because the Trait2Trait module only relies on the predicted traits as inputs, the inferred neural network is small (in terms of numbers of layers, neurons, and parameters). In our application, the module was trained efficiently, taking only 1.5 min for simulated data (compared with 3 min for DL) and 19 min for real data, all on Python using one CPU. Although in the present study the module was implemented and trained using Python, it can be easily implemented in other frameworks (e.g., in R using the Keras library; Kalinowski et al., 2024) and trained on a CPU within minutes. The reason DL-

GBLUP outperformed the DL model might be that GB-LUP is known for accurately predicting additive genetic values, ensuring the initial predictions ($PGV_1$) capture the additive effect with minimal errors. In contrast, the DL SNP2Trait model may face challenges such as overfitting or insufficient capture of additive genetic variance, which could affect its predictive performance. Starting with a precise and comprehensive $PGV_1$ from genomic data ensures accurate modeling and adjustment according to the full relationship between traits using the Trait2Trait module.

Although research on the applications of DL for GP has become more widespread in recent years, works that apply DL for multitrait prediction models have been less popular. Nevertheless, a previous study considered that relationships between the elements of the output layer could be learned and captured automatically by a neural network with shared neurons and weights (Waldmann et al., 2020). However, this study did not explore whether these trait relationships are linear or nonlinear. In contrast, our study reveals that the SNP2Trait module for the
DL model was unable to detect nonlinear relationships. Instead, the DL model required a proper and exclusive network to map one trait (input) to another (output) to capture nonlinearity. A convolutional neural network was tested by including a convolutional block before the MLP block in the SNPs2Trait module. However, this architecture underperformed compared with the MLP model (Supplemental Table S3, see Notes).

In our study on simulated data, the dependent traits were partially conditioned on the GV of a reference trait directly simulated from genomic data. Although the relationships considered in our study are over-characterized compared with real-world data, this controlled setup was essential to develop and validate our DLGBLUP methodology. By testing the model under exaggerated nonlinear relationships, we were able to thoroughly assess its robustness in detecting and capturing such effects.

To further evaluate the model's performance in a more realistic context, we also considered a scenario where the input genomic data included not only QTLs, but also additional noncausative SNPs. In this case, the model remained capable of capturing nonlinear relationships and, in some instances, even improved prediction accuracy. These results demonstrate the method's ability to handle more complex genomic architectures beyond the idealized setting.

Additionally, we employed an ideal prediction scenario in which QTLs were known, leading to very high prediction accuracy—levels that are typically unattainable in real data. This is expected because in simulations, accuracy is measured as the correlation between the predicted GV and TGV, whereas in real data, accuracy is typically computed as the correlation between predicted values and phenotypes, where the GV is not directly observed. As discussed in previous studies (e.g., Frouin et al., 2020), the accuracy in real data is inherently constrained by a function of the trait heritability.

Although we recognize that the purely QTL-based simulation does not fully capture real-world complexities, it provides a necessary proof of concept before applying the model to empirical data. Future work will focus on testing the methodology using real datasets where genetic architectures are less structured and QTLs are unknown.

Regarding the differences between functions, we observed that the logistic function resulted in similar accuracy to GBLUP, whereas the quadratic function led to higher accuracy than GBLUP. This is due to the complexity of the genetic relationships modeled in each case, with quadratic functions capturing stronger deviations from linearity compared with logistic functions, which were close to a linear form. Additionally, in scenarios where causal loci were removed, the superiority of DL-GBLUP was not consistently observed, highlighting the dependency of the method's advantage on the availability of informative markers and the effect of the QTL-to-SNP ratio.

In such cases, the prediction performance of GBLUP and even DL SNP2Trait was considerably reduced, indicating missing genetic information. Because these predicted values serve as inputs to the Trait2Trait module, the model's ability to fully recover and predict nonlinear relationships was also limited.

In dairy cattle, it may be possible to find nonlinear genetic relationships between various traits, although these relationships have not been widely explored, given that usually BLUP is used as the evaluation model, limiting all estimated between-trait relationships to be linear. For example, growth rate and feed efficiency might positively correlate up to a point, after which efficiency plateaus or declines. Fertility and milk production could exhibit a nonlinear relationship, with moderate production levels having minimal impact on fertility, but very high levels reducing it. Longevity and production traits may positively correlate at moderate levels, whereas extreme production reduces lifespan. Mastitis and milking speed may also exhibit nonlinear relationships, with both high and low speeds increasing mastitis risk due to tissue damage or prolonged exposure, respectively. These possibilities suggest that advanced methods may be needed to fully understand and model these complex relationships for optimal breeding and management strategies.

On the real data from the French Holstein dairy cattle, our method was capable of detecting nonlinear relationships between multiple traits (Figure 8). This finding is important for selection, because as shown with simulations, greater genetic progress can be achieved when nonlinear relationships are accounted for, even if the accuracy of the predicted breeding value is not notably changed. A good understanding of these relationships can be useful to manipulate and possibly reshape the relationship through selection in the next generation and improve genetic gain. However, this detection did not improve the prediction accuracy of these traits, which can be justified by the fact that these nonlinear relationships were limited at the population level. Additionally, the contributions of both linear and nonlinear relationships were smaller than those observed in the simulated data. Note that in our evaluation, we used YD obtained from a previous preprocessing as the target, which can alter the actual form of the relationship between the traits and lead to these results. Moreover, the application of such models to raw performance data would need to consider additional factors, such as phenotypic records that may deviate from normality and multiple records for many production traits. The complexity of these models would exceed the scope of our study's objectives, which were to test whether DL could identify nonlinear relationships between traits and use this information to improve prediction accuracy. Further work is required to extend our approach to raw real data.
Better understanding how multiple traits involved in a breeding program are related is a potential way to improve the genetic progress obtained by artificial selection. Here, we maintained the use of a linear SI to select individuals, but considered the additive GV obtained with the different models (GBLUP, DL, and DLGBLUP) for comparison and showed that for all traits after 7 generations, genetic gains from SI based on the additive GV from either DL or DLGBLUP were equal or superior to the gains from SI based on the additive GV from GBLUP. We also found that the form of relationships between traits could be detected in the subsequent generation, indicating that these relationships are inherited and influenced by genetic factors. However, as selection continued over several generations, the relationships, especially for quadratic traits, tended to become more linear. This shift occurs because selection tends to favor a particular subset of the population—in this case, individuals exhibiting linear positive relationships between traits. Over time, selection reduces genetic diversity in the population, leading to a more linear relationship between traits. This phenomenon could explain why our method did not yield greater prediction accuracy in real data. Nonetheless, our proposed method could be useful for evaluating new traits or reassessing previously selected traits when breeding objectives change. In such cases, existing nonlinear effects may reemerge.

## CONCLUSIONS

In this work, we studied the possibility of nonlinear genetic relationships between traits and proposed the hybrid DLGBLUP model to capture them. By integrating GBLUP predictions with a DL module, DLGBLUP can adjust predictions for potential nonlinear relationships when they exist. Applied to simulated data, DLGBLUP improved the accuracy of GV predictions compared to GBLUP in scenarios where nonlinear relationships were present. Moreover, selection based on a SI built from the additive GV predicted by DLGBLUP or DL led to greater genetic gain than using GBLUP alone. Applied to real data from French Holstein cattle, DLGBLUP identified nonlinear genetic relationships between commercial traits, even though prediction accuracy did not improve. These findings suggest that nonlinear relationships between traits can influence selection outcomes and should be considered, particularly when breeding objectives evolve or new traits are introduced.

## NOTES

This research was funded by a CIFRE PhD grant from Eliance (Paris, France), with financial support from the Association Nationale de la Recherche et de la Tech-

nologie (ANRT-Cifre; Paris, France) and APIS-GENE (Paris, France). This work was supported by INRAE Metaprogramme DIGIT-BIO (Digital biology to explore and predict living organisms in their environment; Paris, France). We thank the GenIALearn team for their valuable discussion and contribution, and Jocelyn de Goer De Herve from INRAE/UMR EPIA (Theix, France) for managing the server on which DL work was done. Supplementary material for this article is available at https://doi.org/10.17605/OSF.IO/PS7F2. All source code in this study and simulated data are freely available at https://github.com/fshokor/DLGBLUP. The real data are not publicly available due to privacy. No human or animal subjects were used, so this analysis did not require approval by an Institutional Animal Care and Use Committee or Institutional Review Board. The authors have not stated any conflicts of interest.

**Nonstandard abbreviations used:** BV = breeding value; DL = deep learning; DLGBLUP = hybrid model that combines both DL and GBLUP; GP = genomic prediction; GV = genetic value; LD = linkage disequilibrium; MAE = mean absolute error; MLP = multilayer perceptron; MSE = mean squared error; PGV = predicted genetic value; SI = selection index; ssGBLUP = single-step GBLUP; TGV = true genetic value; YD = yield deviation.

## REFERENCES

* Abdollahi-Arpanahi, R., D. Gianola, and F. Peñagaricano. 2020. Deep learning versus parametric and ensemble methods for genomic prediction of complex phenotypes. *Genet. Sel. Evol.* 52:12. https://doi.org/10.1186/s12711-020-00531-z.

* Akingbuwa, W., and M. Nivard. 2023. Modelling non-linear genetic correlations between BMI and depression. *Eur. Neuropsychopharmacol.* 75(Suppl. 1):S66. https://doi.org/10.1016/j.euroneuro.2023.08.129.

* Alzubaidi, L., J. Zhang, A. J. Humaidi, A. Al-Dujaili, Y. Duan, O. Al-Shamma, J. Santamaría, M. A. Fadhel, M. Al-Amidie, and L. Farhan. 2021. Review of deep learning: concepts, CNN architectures, challenges, applications, future directions. *J. Big Data* 8:53. https://doi.org/10.1186/s40537-021-00444-8.

* Bellot, P., G. de los Campos, and M. Pérez-Enciso. 2018. Can deep learning improve genomic prediction of complex human traits? *Genetics* 210:809–819. https://doi.org/10.1534/genetics.118.301298.

* Calus, M. P., and R. F. Veerkamp. 2011. Accuracy of multi-trait genomic selection using different methods. *Genet. Sel. Evol.* 43:26. https://doi.org/10.1186/1297-9686-43-26.

* Castañeda-Bustos, V. J., H. H. Montaldo, M. Valencia-Posadas, L. Shepard, S. Pérez-Elizalde, O. Hernández-Mendo, and G. Torres-Hernández. 2017. Linear and nonlinear genetic relationships between type traits and productive life in US dairy goats. *J. Dairy Sci.* 100:1232–1245. https://doi.org/10.3168/jds.2016-11313.

* Chen, C., S. A. Bhuiyan, E. Ross, O. Powell, E. Dinglasan, X. Wei, F. Atkin, E. Deomano, and B. Hayes. 2024. Genomic prediction for sugarcane diseases including hybrid Bayesian-machine learning approaches. *Front. Plant Sci.* 15:1398903. https://doi.org/10.3389/fpls.2024.1398903.

* Choudhary, K., B. DeCost, C. Chen, A. Jain, F. Tavazza, R. Cohn, C. W. Park, A. Choudhary, A. Agrawal, S. J. L. Billinge, E. Holm, S. P. Ong, and C. Wolverton. 2022. Recent advances and applications of deep learning methods in materials science. *NPJ Comput. Mater.* 8:59. https://doi.org/10.1038/s41524-022-00734-6.
Costa, W. G. D., M. O. Celeri, I. P. Barbosa, G. N. Silva, C. F. Azevedo, A. Borem, M. Nascimento, and C. D. Cruz. 2022. Genomic prediction through machine learning and neural networks for traits with epistasis. Comput. Struct. Biotechnol. J. 20:5490–5499. https://doi.org/10.1016/j.csbj.2022.09.029.

Cuyabano, B. 2021. *GenEval. Accessed Sep. 4, 2023.* https://github.com/bcuyabano/GenEval.

de Roos, A. P. W., B. J. Hayes, R. J. Spelman, and M. E. Goddard. 2008. Linkage disequilibrium and persistence of phase in Holstein-Friesian, Jersey and Angus cattle. Genetics 179:1503–1512. https://doi.org/10.1534/genetics.107.084301.

Erbe, M., B. J. Hayes, L. K. Matukumalli, S. Goswami, P. J. Bowman, C. M. Reich, B. A. Mason, and M. E. Goddard. 2012. Improving accuracy of genomic predictions within and between dairy cattle breeds with imputed high-density single nucleotide polymorphism panels. J. Dairy Sci. 95:4114–4129. https://doi.org/10.3168/jds.2011-5019.

Frouin, A., C. Dandine-Roulland, M. Pierre-Jean, J.-F. Deleuze, C. Ambroise, and E. Le Floch. 2020. Exploring the link between additive heritability and prediction accuracy from a ridge regression perspective. Front. Genet. 11:581594. https://doi.org/10.3389/fgene.2020.581594.

Fuerst-Waltl, B., A. Essl, and J. Sölkner. 1997. Nonlinear genetic relationships between traits and their implications on the estimation of genetic parameters. J. Anim. Sci. 75:3119–3125. https://doi.org/10.2527/1997.75123119x.

Fuerst-Waltl, B., J. Sölkner, A. Essl, I. Hoeschele, and C. Fuerst. 1998. Non-linearity in the genetic relationship between milk yield and type traits in Holstein cattle. Livest. Prod. Sci. 57:41–47. https://doi.org/10.1016/S0301-6226(98)00161-4.

GenEval. 2024. Évaluation génétique des bovins laitiers – Production laitière, morphologie, caractères fonctionnels. Accessed Apr. 7, 2025. https://www.geneval.fr/_files/ugd/a0004f_640ae3ffc1d446ff80d5406acc14837e.pdf.

Goodfellow, I., Y. Bengio, and A. Courville. 2016. Deep Learning. MIT Press, Cambridge, MA.

Hastie, T., R. Tibshirani, and J. Friedman. 2001. The Elements of Statistical Learning. Springer, New York, NY.

Henderson, C. R., and R. L. Quaas. 1976. Multiple trait evaluation using relatives' records. J. Anim. Sci. 43:1188–1197. https://doi.org/10.2527/jas1976.4361188x.

Hill, W. G. 2013. Genetic correlation. Pages 237–239 in Brenner's Encyclopedia of Genetics. Elsevier, San Diego, CA.

Jubair, S., J. R. Tucker, N. Henderson, C. W. Hiebert, A. Badea, M. Domaratzki, and W. G. D. Fernando. 2021. GPTransformer: A transformer-based deep learning method for predicting Fusarium-related traits in barley. Front. Plant Sci. 12:761402. https://doi.org/10.3389/fpls.2021.761402.

Kalinowski, T., D. Falbel, J. J. Allaire, and F. Chollet. RStudio, Google, Y. Tang, W. Van Der Bijl, M. Studer, and S. Keydana. 2024. keras: R Interface to “Keras.” R package version 2.15.0. Accessed Apr. 7, 2025. https://CRAN.R-project.org/package=keras.

LeCun, Y., Y. Bengio, and G. Hinton. 2015. Deep learning. Nature 521:436–444. https://doi.org/10.1038/nature14539.

Lee, H. J., J. H. Lee, C. Gondro, Y. J. Koh, and S. H. Lee. 2023. deepGBLUP: Joint deep learning networks and GBLUP framework for accurate genomic prediction of complex traits in Korean native cattle. Genet. Sel. Evol. 55:6. https://doi.org/10.1186/s12711-023-00825-y.

Ma, W., Z. Qiu, J. Song, J. Li, Q. Cheng, J. Zhai, and C. Ma. 2018. A deep convolutional neural network approach for predicting phe notypes from genotypes. Planta 248:1307–1318. https://doi.org/10.1007/s00425-018-2976-9.

Meuwissen, T. H., B. J. Hayes, and M. E. Goddard. 2001. Prediction of total genetic value using genome-wide dense marker maps. Genetics 157:1819–1829. https://doi.org/10.1093/genetics/157.4.1819.

Misztal, I., A. Legarra, and I. Aguilar. 2009. Computing procedures for genetic evaluation including phenotypic, full pedigree, and genomic information. J. Dairy Sci. 92:4648–4655. https://doi.org/10.3168/jds.2009-2064.

Montesinos-López, O. A., A. Montesinos-López, J. Crossa, D. Gianola, C. M. Hernández-Suárez, and J. Martín-Vallejo. 2018. Multi-trait,

multi-environment deep learning modeling for genomic-enabled prediction of plant traits. G3 (Bethesda) 8:3829–3840. https://doi.org/10.1534/g3.118.200728.

Montesinos-López, O. A., A. Montesinos-López, J. Crossa, F. H. Toledo, O. Pérez-Hernández, K. M. Eskridge, and J. Rutkoski. 2016. A genomic Bayesian multi-trait and multi-environment model. G3 (Bethesda) 6:2725–2744. https://doi.org/10.1534/g3.116.032359.

Montesinos-López, O. A., A. Montesinos-López, P. Pérez-Rodríguez, J. A. A. Barrón-López, J. W. R. Martini, S. B. Fajardo-Flores, L. S. Gaytan-Lugo, P. C. Santana-Mancilla, and J. Crossa. 2021. A review of deep learning applications for genomic selection. BMC Genomics 22:19. https://doi.org/10.1186/s12864-020-07319-x.

Nazzicari, N., and F. Biscarini. 2022. Stacked kinship CNN vs. GBLUP for genomic predictions of additive and complex continuous phe notypes. Sci. Rep. 12:19889. https://doi.org/10.1038/s41598-022-24405-0.

Pérez-Enciso, M., and L. M. Zingaretti. 2019. A guide for using deep learning for complex trait genomic prediction. Genes (Basel) 10:553. https://doi.org/10.3390/genes10070553.

Pérez-Rodríguez, P., and G. de los Campos. 2022. Multitrait Bayesian shrinkage and variable selection models with the BGLR-R package. Genetics 222:iyac112. https://doi.org/10.1093/genesitics/iyac112.

Pook, T., J. Freudenthal, A. Korte, and H. Simianer. 2020. Using local convolutional neural networks for genomic prediction. Front. Genet. 11:561497. https://doi.org/10.3389/fgene.2020.561497.

Shinde, P., and S. Shah. 2018. A review of machine learning and deep learning applications. Pages 1–6 in Proc. 2018 Fourth International Conference on Computing, Communication, Control and Automation (ICCUBEA), Pune, India. IEEE.

Sölkner, J., and B. Fuerst-Waltl. 1996. Nonlinear heritabilities and curvilinear relationships: Genetic factors complicating selection for functional traits. Interbull Bull. 12:226–233.

Tribout, T., V. Ducrocq, and D. Boichard. 2020. HSSGBLUP: A single-step SNP-BLUP genomic evaluation software adapted to large livestock populations. Proc. International Conference on Quantitative Genetics, Melbourne, Australia.

van Hilten, A., S. A. Kushner, M. Kayser, M. A. Ikram, H. H. H. Adams, C. C. W. Klaver, W. J. Niessen, and G. V. Roshchupkin. 2021. Gen-Net framework: Interpretable deep learning for predicting phenotypes from genetic data. Commun. Biol. 4:1094. https://doi.org/10.1038/s42003-021-02622-z.

van Rheenen, W., W. J. Peyrot, A. J. Schork, S. H. Lee, and N. R. Wray. 2019. Genetic correlations of polygenic disease traits: From theory to practice. Nat. Rev. Genet. 20:567–581. https://doi.org/10.1038/s41576-019-0137-z.

VanRaden, P. M. 2008. Efficient methods to compute genomic predictions. J. Dairy Sci. 91:4414–4423. https://doi.org/10.3168/jds.2007-0980.

Vaswani, A., N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin. 2017. Attention is all you need. Adv. Neural Inf. Process. Syst. 30:5998–6008.

Waldmann, P., C. Pfeiffer, and G. Mészáros. 2020. Sparse convolutional neural networks for genome-wide prediction. Front. Genet. 11:25. https://doi.org/10.3389/fgene.2020.00025.

Wang, K., M. A. Abid, A. Rasheed, J. Crossa, S. Hearne, and H. Li. 2023. DNNGP, a deep neural network-based method for genomic prediction using multi-omics data in plants. Mol. Plant 16:279–293. https://doi.org/10.1016/j.molp.2022.11.004.

Zingaretti, L. M., S. A. Gezan, L. F. Ferrão, L. F. Osorio, A. Monfort, P. R. Muñoz, V. M. Whitaker, and M. Pérez-Enciso. 2020. Exploring deep learning for complex trait genomic prediction in polyploid outcrossing species. Front. Plant Sci. 11:25. https://doi.org/10.3389/fpls.2020.00025.

## ORCIDS

F. Shokor, https://orcid.org/0009-0007-5315-1530

T. Tribout, https://orcid.org/0000-0003-2500-7375

B. C. D. Cuyabano, https://orcid.org/0000-0003-0060-8879