ADAM-multi: software to simulate complex breeding programs for animals and plants with different ploidy levels and generalized genotypic effect models to account for multiple alleles

Stochastic simulation software, ADAM, has been developed for the purpose of breeding optimization in animals and plants, and for validation of statistical models used in genetic evaluations. Just like other common simulation programs, ADAM assumed the bi-allelic state of quantitative trait locus (QTL). While the bi-allelic state of marker loci is due to the common choice of genotyping technology of single nucleotide polymorphism (SNP) chip, the assumption may not hold for the linked QTL. In the version of ADAM-Multi, we employ a novel simulation model capable of simulating additive, dominance, and epistatic genotypic effects for species with different levels of ploidy, providing with a more realistic assumption of multiple allelism for QTL variants. When assuming bi-allelic QTL, our proposed model becomes identical to the model assumption in common simulation programs, and in genetic textbooks. Along with the description of the updated simulation model in ADAM-Multi, this paper shows two small-scale studies that investigate the effects of multi-allelic versus bi-allelic assumptions in simulation and the use of different prediction models in a single-population breeding program for potatoes. We found that genomic models using dense bi-allelic markers could effectively predicted breeding values of individuals in a well-structure population despite the presence of multi-allelic QTL. Additionally, the small-scale study indicated that including non-additive genetic effects in the prediction model for selection did not lead to an improvement in the rate of genetic gains of the breeding program.

1 Introduction

Stochastic simulation is a cost-effective and powerful tool to optimize breeding programs with reduced experimental costs. Such a tool unlocks possibilities for investigating alternative breeding schemes, in order to maximizing genetic gains of the breeding program at a given input of resources. Software package, ADAM (Pedersen et al., 2009; Liu et al., 2019), has been developed for the purpose of breeding optimization in pig, cattle, fish and plants (Bengtsson et al., 2022; Tessema et al., 2020; Zaalberg et al., 2022; Chu et al., 2020). The tool is also very useful in validating statistical genetic prediction models (Romé et al., 2023; Chu et al., 2021) and in studying methods for preserving genetic diversity in breeding programs (Henryon et al., 2015). Over time the software has been further developed and updated with many features including extension to non-additive genetic models (Chu et al., 2024), indirect genetic effects (Chu et al., 2021), categorical traits (Gebreyesus et al., 2020), definition of true inbreeding (Henryon et al., 2019), or extension to plant breeding (Liu et al., 2019).

ADAM (Pedersen et al., 2009) simulates genotypic effects for each allele in a quantitative trait locus (QTL). This model is unique from other software like AlphaSim (Gaynor et al., 2021; Faux et al., 2016), ChromaX (Younis et al., 2023), MoBPS (Pook et al., 2020), SeqBreed (Pérez-Enciso et al., 2020) and XSim (Chen et al., 2022) that use substitution genotypic effects of QTL. However, all of these tools including ADAM assume a bi-allelic state of each segregating QTL. This assumption might come from current, common genetic models (Christensen et al., 2012; VanRaden, 2008; Falconer and Mackay, 1996) that assume substitution effects of one allele to its alternative. While the bi-allelic state of marker loci is due to the common choice of genotyping technology that yield single nucleotide polymorphism (SNP). The assumption of bi-allelic QTL may not true for all QTL. Multiple alleles have been shown in numerous QTL (Biová et al., 2024; Jiang et al., 2020). In addition, multi-allelic models of QTL are more reasonable explanations for different functional genetic effects from unrelated populations (González-Diéguez et al., 2021).

Multi-allelic models have been developed for genomic prediction in diploid species (Álvarez-Castro and Crujeiras, 2019; Álvarez-Castro and Yang, 2011; Yang and Álvarez-Castro, 2008; Da, 2015). Relevant model based on haplotype blocks also have been shown for genomic prediction (Weber et al., 2023). Thérèse Navarro et al. (2022) has developed a package for genome wide association studies (GWAS) of polyploid populations with multi-allelic models, but only the additive genetic effects were included in the model. For the purpose of simulation, however, we are not aware of any studies that have used multi-allelic models with additive, dominance and epistatic genetic effects, or accounting for different levels of ploidy.

In addition to bi-allelic assumptions, earlier version of ADAM (Pedersen et al., 2009; Liu et al., 2019) use was limited to diploid species only. Many economically important species like potato, banana, sugar cane and some fish orders of salmonids and common carps are polyploid. Extension of simulation models to different ploidy levels is necessary for ADAM (Pedersen et al., 2009; Liu et al., 2019) to design complex breeding schemes for these species.

This paper will describe new features of software package, now called ADAM-Multi, for simulating breeding programs for plants and animals. The focus will be on description of genotypic models for simulating traits with genotypic effects of additive, dominance, and epistatic genetics for species with different ploidy levels including extensions to multi-allelic assumptions. The methods implemented are illustrated in two examples that study the effects of different assumptions on number of alleles, ploidy level and different prediction models used in selection.

2 Materials and methods

2.1 Genotypic models for simulation

We aim to simulate genotypic effects that are generalized to multi-allelic QTL with number ($n_B$) of alleles, and the genome with ploidy level of $n_{ploidy}$. Assuming a QTL with alleles $B_1$, $B_2$,…, $B_{n_B}$, ADAM-Multi uses following model to simulate additive genotypic value ($a$) at one locus:

$a={\displaystyle \sum _{i_B}^{n_B}}{t}_{i_B}^a{a}_{i_B}=t_1^a{a}_1+t_2^a{a}_2+\dots +t_{n_B}^a{a}_{n_B}(1)$

where $a$ is the additive genotypic value of a QTL; $a_{i_B}$ is genotypic additive effect of allele $B_{i_B}$ (or called $i_B$ for short) at the QTL; $t_{i_B}^a$ is the additive covariate for allele $i_B$ that is scaled genotype dosage calculated as in the AlphaSimR software (Gaynor et al., 2021):

$t_{i_B}^a=\left(t_{i_B}-\frac{n_{ploidy}}{2}\right)\left(\frac{2}{n_{ploidy}}\right)(2)$

where $t_{i_B}$ is a raw genotype dosage, or number of copies of allele $i_B$ at the locus, $n_{ploidy}$ is the ploidy level of genome. Key notations are defined in Table 1.

Table 1

Table 1. List of key symbols.

Similarly, the genotypic model for dominance value ($d$) of a QTL at the locus level is:

$d={\displaystyle \sum _{i_B}^{n_B}}{t}_{i_B}^d{d}_{i_B}=t_1^d{d}_1+t_2^d{d}_2+\dots +t_{n_B}^d{d}_{n_B}(3)$

where $d_{i_B}$ is dominance genotypic effect of allele $i_B;t_{i_B}^d$ is the dominance covariate for allele $i_B$ that is scaled genotype dosage calculated as in AlphaSimR (Gaynor et al., 2021):

$t_{i_B}^d=t_{i_B}\left(n_{ploidy}-t_{i_B}\right){\left(\frac{2}{n_{ploidy}}\right)}^2(4)$

This simulation model assumes digenic dominance for each allele, i.e., each allele has the same dominance effect with all other alleles. Table 2 shows examples of $t$, $t^a$ and $t^d$ for diploid and tetraploid genome assuming $n_B=2$ with allele $B_1$ and $B_2$.

Table 2

Table 2. Conversion from raw genotype dosages (${\mathit{t}}_{\mathbf{1}}$, ${\mathit{t}}_{\mathbf{2}}$) to additive (${\mathit{t}}_{\mathbf{1}}^{\mathit{a}}$, ${\mathit{t}}_{\mathbf{2}}^{\mathit{a}}$) and dominance (${\mathit{t}}_{\mathbf{1}}^{\mathit{d}}$, ${\mathit{t}}_{\mathbf{2}}^{\mathit{d}}$) covariates when assuming bi-allelic loci.

The simulation model for additive × additive genotypic value ${\left(aa\right)}^{kl}$ of the two-locus epistatic interaction between the pair of loci $k$ and $l$:

${\left(aa\right)}^{kl}=\left\{\left[\begin{array}{cccc}{t}_1^{a,l}& t_2^{a,l}& \dots & t_{n_B}^{a,l}\end{array}\right]\otimes \left[\begin{array}{cccc}{t}_1^{a,k}& t_2^{a,k}& \dots & t_{n_B}^{a,k}\end{array}\right]\right\}\cdot \left[\begin{array}{c}{\left(aa\right)}_1^{kl}\\ {\left(aa\right)}_2^{kl}\\ \dots \\ {\left(aa\right)}_{n_B\times n_B}^{kl}\end{array}\right](5)$

where $\left[\begin{array}{cccc}{t}_1^{a,x}& t_2^{a,x}& \dots & t_{n_B}^{a,x}\end{array}\right]$ is a vector of additive covariates for locus $x$ ($k$ or $l$) with $n_B$ elements; $\otimes$ denotes the Kronecker product; denotes the symbol for matrix multiplication; $\left[\begin{array}{c}{\left(aa\right)}_1^{kl}\\ {\left(aa\right)}_2^{kl}\\ \dots \\ {\left(aa\right)}_{n_B\times n_B}^{kl}\end{array}\right]$ is vector of additive × additive genotypic effects that have $n_B\times n_B$ elements. In Equation 5, the number of alleles in loci $k$ and $l$ are the same, and equal to $n_B$. Simulation models in this paper consider a fixed number $n_B$ for all QTL, even if not all alleles in a QTL are segregating. Genotypic models in case of arbitrary number of alleles for QTL that set the effects of non-segregating alleles to zero can be found in Supplementary Appendix 1.

The simulation model for additive-dominance genotypic value ${\left(ad\right)}^{kl}$ of the epistatic interaction is the sum of additive × dominance and dominance × additive interaction between the pair of loci $k$ and $l$:

${\left(ad\right)}^{kl}=\left\{\left[\begin{array}{cccc}{t}_1^{d,l}& t_2^{d,l}& \dots & t_{n_B}^{d,l}\end{array}\right]\otimes \left[\begin{array}{cccc}{t}_1^{a,k}& t_2^{a,k}& \dots & t_{n_B}^{a,k}\end{array}\right]\right\}\cdot \left[\begin{array}{c}{\left(ad\right)}_1^{kl}\\ {\left(ad\right)}_2^{kl}\\ \dots \\ {\left(ad\right)}_{n_B\times n_B}^{kl}\end{array}\right]+\left\{\left[\begin{array}{cccc}{t}_1^{a,l}& t_2^{a,l}& \dots & t_{n_B}^{a,l}\end{array}\right]\otimes \left[\begin{array}{cccc}{t}_1^{d,k}& t_2^{d,k}& \dots & t_{n_B}^{d,k}\end{array}\right]\right\}\cdot \left[\begin{array}{c}{\left(da\right)}_1^{kl}\\ {\left(da\right)}_2^{kl}\\ \dots \\ {\left(da\right)}_{n_B\times n_B}^{kl}\end{array}\right](6)$

where $\left[\begin{array}{cccc}{t}_1^{d,x}& t_2^{d,x}& \dots & t_{n_B}^{d,x}\end{array}\right]$ is a vector of dominance covariates (Equation 4) for locus $x$ with $n_B$ elements; $\left[\begin{array}{c}{\left(ad\right)}_1^{kl}\\ {\left(ad\right)}_2^{kl}\\ \dots \\ {\left(ad\right)}_{n_B\times n_B}^{kl}\end{array}\right]\text{\hspace{0.17em}and\hspace{0.17em}}\left[\begin{array}{c}{\left(da\right)}_1^{kl}\\ {\left(da\right)}_2^{kl}\\ \dots \\ {\left(da\right)}_{n_B\times n_B}^{kl}\end{array}\right]$ are vectors of additive × dominance (${\left(ad\right)}_x^{kl}$), and dominance × additive (${\left(da\right)}_x^{kl}$) genotypic effects that have $n_B\times n_B$ elements. The value ${\left(ad\right)}_x^{kl}$ is different from ${\left(da\right)}_x^{kl}$.

The simulation model for dominance × dominance genotypic value ${\left(dd\right)}^{kl}$ of the epistatic interaction between the pair of loci $k$ and $l$:

${\left(dd\right)}^{kl}=\left\{\left[\begin{array}{cccc}{t}_1^{d,l}& t_2^{d,l}& \dots & t_{n_B}^{d,l}\end{array}\right]\otimes \left[\begin{array}{cccc}{t}_1^{d,k}& t_2^{d,k}& \dots & t_{n_B}^{d,k}\end{array}\right]\right\}\cdot \left[\begin{array}{c}{\left(dd\right)}_1^{kl}\\ {\left(dd\right)}_2^{kl}\\ \dots \\ {\left(dd\right)}_{n_B\times n_B}^{kl}\end{array}\right](7)$

where $\left[\begin{array}{c}{\left(dd\right)}_1^{kl}\\ {\left(dd\right)}_2^{kl}\\ \dots \\ {\left(dd\right)}_{n_B\times n_B}^{kl}\end{array}\right]$ is vector of dominance × dominance genotypic effects that have $n_B\times n_B$ elements.

So far, the simulated genotypic values were presented at the levels of locus and loci pairs. Here, we present the model for simulating genotypic value at the individual level, which is the sum effects of all QTLs and epistatic pair interactions. The model in a matrix form for total genotypic value $g_i$ of individual $i$ is:

$g_i=\text{tr}\left({\mathbf{t}}_{\mathbf{i}}^{\mathbf{a}}\cdot \mathbf{a}\right)+\text{tr}\left({\mathbf{t}}_{\mathbf{i}}^{\mathbf{d}}\cdot \mathbf{d}\right)+\text{tr}\left[{\mathbf{t}}_{\mathbf{i}}^{\mathbf{\text{aa}}}\cdot \left(\mathbf{\text{aa}}\right)\right]+\text{tr}\left[{\mathbf{t}}_{\mathbf{i}}^{\mathbf{\text{dd}}}\cdot \left(\mathbf{\text{dd}}\right)\right]+\text{tr}\left[{\mathbf{t}}_{\mathbf{i}}^{\mathbf{\text{ad}}}\cdot \left(\mathbf{\text{ad}}\right)\right]+\text{tr}\left[{\mathbf{t}}_{\mathbf{i}}^{\mathbf{\text{da}}}\cdot \left(\mathbf{\text{da}}\right)\right](8)$

where $\mathbf{a}$ is a $n_B\times n_{qtl}$ matrix:

$\mathbf{a}=\left[\begin{array}{cccc}{a}_1^{j_1}& a_1^{j_2}& \dots & a_1^{n_{qtl}}\\ a_2^{j_1}& a_2^{j_2}& \dots & a_2^{n_{qtl}}\\ \dots & \dots & \dots & \dots \\ a_{n_B}^{j_1}& a_{n_B}^{j_2}& \dots & a_{n_B}^{n_{qtl}}\end{array}\right],$

where $n_{qtl}$ is the number of $\text{QTL};a_{i_B}^{j_{qtl}}$ is the additive genotypic effect of allele $i_B$ at $\text{QTL\hspace{0.17em}}{j}_{qtl};{\mathbf{t}}_{\mathbf{i}}^{\mathbf{a}}\text{\hspace{0.17em}is\hspace{0.17em}}\mathrm{a}{n}_{qtl}\times n_B$ matrix:

${\mathbf{t}}_{\mathbf{i}}^{\mathbf{a}}=\left[\begin{array}{cccc}{t}_{1,i}^{a,1}& t_{2,i}^{a,1}& \dots & t_{n_B,i}^{a,1}\\ t_{1,i}^{a,2}& t_{2,i}^{a,2}& \dots & t_{n_B,i}^{a,2}\\ \dots & \dots & \dots & \dots \\ t_{1,i}^{a,n_{qtl}}& t_{2,i}^{a,n_{qtl}}& \dots & t_{n_B,i}^{a,n_{qtl}}\end{array}\right],$

where $t_{i_B,i}^{a,j_{qtl}}$ is the additive covariate of allele $i_B$ at $\text{QTL\hspace{0.17em}}{j}_{qtl}$ of individual $i.t_{i_B,i}^{a,j_{qtl}}$ can be calculated based on the genotype of individual $i$ at locus $j_{qtl}$ using Equation 2; $\text{tr}\left(\right)$ is the trace of a matrix. Similarly, $\mathbf{d}$ is a $n_B\times n_{qtl}$ matrix of dominance genotypic effects; ${\mathbf{t}}_{\mathbf{i}}^{\mathbf{d}}$ is a $n_{qtl}\times n_B$ matrix of dominance covariate of individual $i$. Matrix $\left(\mathbf{\text{aa}}\right)$ has a dimension of $\left(n_B\times n_B\right)$ rows and $n_{ep}$ columns:

$\left(\mathbf{\text{aa}}\right)=\left[\begin{array}{cccc}{\left(aa\right)}_{1,1}^{kl}& {\left(aa\right)}_{2,1}^{kl}& \dots & {\left(aa\right)}_{n_{ep},1}^{kl}\\ {\left(aa\right)}_{1,2}^{kl}& {\left(aa\right)}_{2,2}^{kl}& \dots & {\left(aa\right)}_{n_{ep},2}^{kl}\\ \dots & \dots & \dots & \dots \\ {\left(aa\right)}_{1,\left(n_B\times n_B\right)}^{kl}& {\left(aa\right)}_{2,\left(n_B\times n_B\right)}^{kl}& \dots & {\left(aa\right)}_{n_{ep},\left(n_B\times n_B\right)}^{kl}\end{array}\right],$

where $n_{ep}$ is number of epistatic interactions between pairs of loci;${\left(aa\right)}_{j_{ep},\left(i_B^k,i_B^l\right)}^{kl}$ is the epistatic additive × additive effect at the interaction pair $j_{ep}$ between allele $i_B^k$ of locus $k$ and allele $i_B^l$ of locus $l.{\mathbf{t}}_{\mathbf{i}}^{\mathbf{\text{aa}}}$ is a matrix with $n_{ep}$ rows and $\left(n_B\times n_B\right)$ columns.Row $j_{ep}$ th of ${\mathbf{t}}_{\mathbf{i}}^{\mathbf{\text{aa}}}$ is set up as in Equation 5. For example,

${\mathbf{t}}_{{\mathbf{j}}_{\mathbf{\text{ep}}},\mathbf{i}}^{\mathbf{\text{aa}}}=\left[\begin{array}{cccc}{t}_{j_{ep},1,i}^{a,l}& t_{j_{ep},2,i}^{a,l}& \dots & t_{j_{ep},n_B,i}^{a,l}\end{array}\right]\otimes \left[\begin{array}{cccc}{t}_{j_{ep},1,i}^{a,k}& t_{j_{ep},2,i}^{a,k}& \dots & t_{j_{ep},n_B,i}^{a,k}\end{array}\right],$

where $t_{j_{ep},i_B^x,i}^{a,x}$ is the additive covariate of allele $i_B^k$ of QTL $k$ or allele $i_B^l$ of QTL $l$ at the interaction pair $j_{ep}$ for individual $i$.

Similarly, $\left(\mathbf{\text{dd}}\right)$ is a $\left(n_B\times n_B\right)\times n_{ep}$ matrix of dominance × dominance genotypic effects; ${\mathbf{t}}_{\mathbf{i}}^{\mathbf{\text{dd}}}\text{\hspace{0.17em}is\hspace{0.17em}}\mathrm{a}{n}_{ep}\times \left(n_B\times n_B\right)$ matrix of dominance × dominance covariates of individual $i$. The genotypic additive-dominance effects consist of two components: additive × dominance $\text{tr}\left[{\mathbf{t}}_{\mathbf{i}}^{\mathbf{\text{ad}}}\cdot \left(\mathbf{\text{ad}}\right)\right]$ and dominance × additive $\text{tr}\left[{\mathbf{t}}_{\mathbf{i}}^{\mathbf{\text{da}}}\cdot \left(\mathbf{\text{da}}\right)\right]$. Similar to $\left(\mathbf{\text{dd}}\right)$, matrices $\left(\mathbf{\text{ad}}\right)$ and $\left(\mathbf{\text{da}}\right)$ have dimension of $\left(n_B\times n_B\right)\times n_{ep}$. ${\mathbf{t}}_{\mathbf{i}}^{\mathbf{\text{ad}}}$ and ${\mathbf{t}}_{\mathbf{i}}^{\mathbf{\text{da}}}$ are matrices with $n_{ep}$ rows and $\left(n_B\times n_B\right)$ columns. Row $j_{ep}$ th of ${\mathbf{t}}_{\mathbf{i}}^{\mathbf{\text{ad}}}$ is set up as in Equation 6. For example,

${\mathbf{t}}_{{\mathbf{j}}_{\mathbf{\text{ep}}},\mathbf{i}}^{\mathbf{\text{ad}}}=\left[\begin{array}{cccc}{t}_{j_{ep},1,i}^{d,l}& t_{j_{ep},2,i}^{d,l}& \dots & t_{j_{ep},n_B,i}^{d,l}\end{array}\right]\otimes \left[\begin{array}{cccc}{t}_{j_{ep},1,i}^{a,k}& t_{j_{ep},2,i}^{a,k}& \dots & t_{j_{ep},n_B,i}^{a,k}\end{array}\right],$

where $t_{j_{ep},i_B^x,i}^{d,x}$ is the dominance covariate of allele $i_B^k$ of QTL $k$ or allele $i_B^l$ of QTL $l$ at the interaction pair $j_{ep}$ for individual $i$. Row $j_{ep}$ th of ${\mathbf{t}}_{\mathbf{i}}^{\mathbf{\text{da}}}$ is:

${\mathbf{t}}_{{\mathbf{j}}_{\mathbf{\text{ep}}},\mathbf{i}}^{\mathbf{\text{da}}}=\left[\begin{array}{cccc}{t}_{j_{ep},1,i}^{a,l}& t_{j_{ep},2,i}^{a,l}& \dots & t_{j_{ep},n_B,i}^{a,l}\end{array}\right]\otimes \left[\begin{array}{cccc}{t}_{j_{ep},1,i}^{d,k}& t_{j_{ep},2,i}^{d,k}& \dots & t_{j_{ep},n_B,i}^{d,k}\end{array}\right].$

2.2 Stochastic simulation

Details of simulation steps and theoretical principles of ADAM-Multi can be found in previous version by Pedersen et al. (2009) and Liu et al. (2019). These principles are also similar to those in AphaSim (Gaynor et al., 2021). Simulation of genomic models with ADAM first starts with founder haplotypes of a defined genome structure. To create linkage disequilibrium (LD) between QTL and markers, ADAM-Multi can be used in case of multi-allelism. Other packages such as QMSIM (Sargolzaei and Schenkel, 2009) and AlphaSim (Gaynor et al., 2021) do not support multi-allelic models, but they can be used to generate the genome with a specified degree of LD in case of bi-allelic loci. The genotypic effects of alleles in QTLs are sampled, and then centered and scaled to user-defined parameters using the founders’ QTL haplotypes (Chu et al., 2024). Steps for generating additive effects of alleles $\mathbf{a}$ ($n_B\times n_{qtl}$ matrix) in ADAM-Multi are:

- Sampling: Each element of matrix ${\mathbf{a}}^{\mathrm{s}}$ is sampled from a user-defined normal distribution, e.g., mean of zero and additive variance ${\sigma }_A^2$ [or N(0, ${\sigma }_A^2$)]. Matrix ${\mathbf{a}}^{\mathrm{s}}$ with the same dimension as $\mathbf{a}$ contains starting values of additive effects.

- Centering: Based on matrix ${\mathbf{a}}^{\mathrm{s}}$ and genotypes of a founder population, we can calculate population mean at each QTL locus. Additive effects of alleles within each QTL (each row of ${\mathbf{a}}^{\mathrm{s}}$) are centered to achieve population mean of zero. For example,: ${a_{i_B}^{j_{qtl}}}^{*}=a_{i_B}^{j_{qt{l}^s}}-{{\mu }^{j_{qtl}}}^s,\text{where\hspace{0.17em}}{{\mu }^{j_{qtl}}}^s$ is the population mean at locus $j_{qtl}$ given ${\mathbf{a}}^{\mathrm{s}}$ (before being centered) and genotypes of the founder population; ${a_{i_B}^{j_{qtl}}}^{*}$ is the prior value after centering. Matrix ${\mathbf{a}}^{*}$ (same dimension as $\mathbf{a}$) with elements of ${a_{i_B}^{j_{qtl}}}^{*}$ is the prior values of additive effects after centering.

- Rescaling: Prior variance ${\sigma }_{A^{*}}^2$ can also be calculated as we know all individuals’ genotype in the population and functional effects of QTL (${\mathbf{a}}^{*}$). Calculation of the variance can be done by different methods including variance by locus, by chromosome, or by individual as in Chu et al. (2024). Additive effects $\mathbf{a}$ are calculated by rescaling prior effects ${\mathbf{a}}^{*}$ to achieve the user-defined variance input ${\sigma }_A^2$ for the founder population. For example:

$\mathbf{a}={\mathbf{a}}^{*}\times \sqrt{\frac{{\sigma }_A^2}{{\sigma }_{A^{*}}^2}}$

The calculated variances of ${\sigma }_A^2$ and ${\sigma }_{A^{*}}^2$ in the rescaling step are functional, biological or genotypic parameters, which are different from classical, statistical quantitative parameters. The differences between functional and statistical variances are detailed and explained in Álvarez-Castro and Carlborg (2007), Álvarez-Castro and Yang (2011), Chu et al. (2024), and Vitezica et al. (2017). Functional effects of dominance and epistasis are also generated based on user-defined inputs of functional variances. Steps for generating dominance effects of alleles $\mathbf{d}$ ($n_B\times n_{qtl}$ matrix) are:

- Sampling: A dominance degree ${\delta }_{i_B}^{d,j_{qtl}}$ for allele $i_B$ of QTL $j_{qtl}$ is sampled from a user-defined normal distribution $N\left({\mu }_{\delta },{\sigma }_{\delta }^2\right)$, where ${\mu }_{\delta }$ is the mean of dominance degree, and ${\sigma }_{\delta }^2$ is the variance of dominance degree. If a user does not provide this distribution, ADAM-Multi default values are $N\left(0.19,0.097\right)$ as in Wellmann and Bennewitz, (2011). Each element of matrix ${\mathbf{d}}^{*}$ (starting (or prior) values of dominance effects) is generated as:

${d_{i_B}^{j_{qtl}}}^{*}={\delta }_{i_B}^{d,j_{qtl}}*\left|a_{i_B}^{j_{qtl}}\right|.$

- Rescaling: This is done similarly to rescaling as for simulating additive effects. Based on matrix ${\mathbf{d}}^{*}$ and QTL genotypes of a founder population, prior variance ${\sigma }_{D^{*}}^2$ can be calculated. Dominance effects $\mathbf{d}$ are the prior effects ${\mathbf{d}}^{*}$ that are rescaled to achieve the user-defined variance inputs of ${\sigma }_D^2$ for the founder population.

Steps for sampling additive × additive effects are.

- Sampling: Each element of matrix ${\left(\mathbf{\text{aa}}\right)}^{*}$ [starting (or prior) values of additive × additive effects] is sampled from a user-defined normal distribution N(0, ${\sigma }_{AA}^2$).

- Rescaling: Prior variance ${\sigma }_{{AA}^{*}}^2$ can be calculated as we know all individuals’ genotype in the population and functional effects of QTL ${\left(\mathbf{\text{aa}}\right)}^{*}$. Calculation of the variance can be done by different methods including variance by pair loci, or by individual as in Chu et al. (2024). Additive × additive effects $\left(\mathbf{\text{aa}}\right)$ are obtained by rescaling prior effects ${\left(\mathbf{\text{aa}}\right)}^{*}$ to achieve the desired variance inputs of ${\sigma }_{AA}^2$ for the founder population.

Steps for sampling dominance × dominance effects are similar to those for additive × additive effects. Steps for sampling additive-dominance effects are:

- Sampling: Each element of matrices ${\left(\mathbf{\text{ad}}\right)}^{*}$ and ${\left(\mathbf{\text{da}}\right)}^{*}$ [starting (or prior) values of additive × dominance and dominance × additive effects] is sampled from a user-defined normal distribution N(0, ${\sigma }_{AD}^2$).

- Rescaling: Prior variance ${\sigma }_{{AD}^{*}}^2$ can be calculated as we know all individuals’ genotype in the population and functional effects of QTL ${\left(\mathbf{\text{ad}}\right)}^{*}$ and ${\left(\mathbf{\text{da}}\right)}^{*}$. Note that the additive-dominance effects of an individuals’ genotype is the sum of additive × dominance and dominance × additive. Effects $\left(\mathbf{\text{ad}}\right)$ and $\left(\mathbf{\text{da}}\right)$ are obtained by rescaling prior effects ${\left(\mathbf{\text{ad}}\right)}^{*}$ and ${\left(\mathbf{\text{da}}\right)}^{*}$ to achieve the desired variance inputs of ${\sigma }_{AD}^2$ for the founder population. For example:

$\left(\mathbf{\text{ad}}\right)={\left(\mathbf{\text{ad}}\right)}^{*}\times \sqrt{\frac{{\sigma }_{AD}^2}{{\sigma }_{{AD}^{*}}^2}}\text{\hspace{0.17em}and\hspace{0.17em}}\left(\mathbf{\text{da}}\right)={\left(\mathbf{\text{da}}\right)}^{*}\times \sqrt{\frac{{\sigma }_{AD}^2}{{\sigma }_{{AD}^{*}}^2}}.$

A centering step is not included for dominance and epistatic effects, but the total genetic value $g_i$ is centered to achieve the population mean of zero based on the founder population. For example, the model for simulating individuals’ phenotypes is:

$y_i=\mu +g_i+e_i(9)$

where $g_i$ is the genetic values of individual $i$ that is constructed as in (Equation 8); $\mu$ is the population mean to re-adjust the mean of the founder population to zero.

In ADAM-Multi, model (Equation 9) can be extended to repeated records, inclusion of non-genetic effects, and multiple traits. The functional genetic effects of additive, dominance, and epistasis are independent of allele frequencies, and in the simulation, they are kept constant across generations. Modeling of genetic recombination during meiosis uses bivalent chromosome pairing (Voorrips and Maliepaard, 2012). A breeding scheme is simulated by combining series of actions: mating, reproduction, phenotyping, genotyping, prediction of breeding values and different selection methods. The use of ADAM-Multi is demonstrated in two examples that study the effects of multi-allelic versus bi-allelic assumptions and the use of different prediction models on accuracy of prediction and genetic gains of breeding programs for potato.

2.3 Example 1

The example uses a simulation model with additive effects only for a single-population breeding scheme. The investigated factors include different multi-allelic assumptions (6 levels) and two levels of ploidy (Table 3). In total, there were 6 × 2 = 12 scenarios simulated in example 1.

Table 3

Table 3. Overview of factors investigated in example 1 and 2.

The simulation model for individuals’ phenotype in example 1 was:

$y_i=\mu +{\displaystyle \sum _{i_B}^{n_B}}{t}_{i_B,i}^{a,1}{a}_{i_B}^{j_1}+{\displaystyle \sum _{i_B}^{n_B}}{t}_{i_B,i}^{a,2}{a}_{i_B}^{j_2}+\dots +{\displaystyle \sum _{i_B}^{n_B}}{t}_{i_B,i}^{a,n_{qtl}}{a}_{i_B}^{n_{qtl}}+e_i(10)$

Simulation of genetic values in (Equation 10) is an extension of (Equation 1) to the sum effects of $n_{qtl}$, or equivalent to the additive part of (Equation 8). Notations and symbols are the same as in (Equation 8). The environment term $e_i$ was drawn from a normal distribution N(0, ${\sigma }_e^2$) with ${\sigma }_e^2=2$. This example used a simplified simulation model as compared to Equation 9 because we would like to assess accuracy of predicted breeding values in selection. When non-additive genetics are included in the model, the definition of accuracy of predicted breeding values with multi-allelic assumption is unclear in literature.

Figure 1 shows the simulation pipeline for this example. The genome of founders was simulated to form LD between QTL and markers using a Fisher-Wright inheritance model (Fisher, 1930). The LD genome consisted of 12 chromosomes with genetic distances emulating that of potato (Massa et al., 2018). The total genome length was 888.6 cM. The initial genome that had 50k marker loci and 10k QTL with an equal frequency for each allele. A historical population with effective population size of 200 was simulated for 1,000 generations of random mating, a simulated bottleneck, and an inheritance pattern of standard Mendelian principles and bivalent chromosome pairing (Voorrips and Maliepaard, 2012). In example 1, 12 founder populations were created corresponding to 6 levels of multi-allelic assumptions and 2 levels of ploidy (Table 3). The 100-individual founders for each of the populations was created, where the genome of these founders consisted of 2k QTL and 10k bi-allelic markers. The QTL and markers were drawn randomly from loci that were segregating with a minor allele frequency ≥0.05. For markers, bi-allelic state was assumed in all scenarios. The percentage of QTL having segerating bi-, tri-, or quad-alleles was corresponding to the assumption of the scenario. In case of multiple allellisms, for example, quad-allelic assumption for a QTL, each of the four alleles must have a minor allele frequency of at least 0.05. The LD pattern in case of bi-allelic loci for diploid and tetraploid genomes can be found in Supplementary Appendix 2. The functional additive variance of the founder populations was simulated at ${\sigma }_A^2=1.0$.

Figure 1

Figure 1. An overview of simulation pipeline in example 1 and 2.

A simplified breeding scheme was simulated for 13 discrete generations. In a generation, 80 parents were crossed, pseudo-randomly with no self-pollination, to create 80 families. A parent could mate any other parents, but each parent could contribute to only maximum of four crosses. Each family had 10 full-sib offspring, thus in total there were 800 offspring per generation. From generation 1 to 5, the 80 parents were randomly selected from the 800 offspring. In generation 6 to 13, the selection of the offspring to be parents in the following generation was based on predicted breeding values. At generation 6, variance components were estimated when the phenotype data consisted of 4,000 individuals. This estimation of variance components ensures that extra variation due to unknown variance components were taken into account. These variance component estimates were used in the models for prediction of breeding values in the subsequent generations. Selection of 80 new parents from 800 individuals were carried out based on the genetic evaluation after the phenotypes in a generation were obtained.

Model Equation 11 were used in example 1 for variance component estimation and prediction of breeding values is as follows:

$\mathbf{y}=\mathbf{\text{Xb}}+{\mathbf{Z}}_{\mathbf{u}}\mathbf{u}+\mathbf{e}(11)$

where $\mathbf{y}$ is the vector of individual phenotypes; $\mathbf{b}$ is a vector of the fixed effects of individuals’ generation; $\mathbf{u}$ is the vector of breeding values $\mathbf{u}\sim \mathbf{N}\left(\mathbf{0},{\mathbf{G}}_{\mathrm{u}}{\sigma }_u^2\right)$, where the relationship matrix for additive genetics (${\mathbf{G}}_{\mathrm{u}}$) was constructed based on bi-allelic markers using method (VanRaden, 2008) for different levels of ploidy assumed in the scenario. The computation of ${\mathbf{G}}_{\mathrm{u}}$ was carried out using AGHmatrix R package (Amadeu et al., 2023). $\mathbf{X}$ and ${\mathbf{Z}}_{\mathrm{u}}$ are design matrices relating individuals to fixed effects and additive genetic effects, respectively. Vector $\mathbf{e}$ is an environmental residual term: $\mathbf{e}\sim \mathbf{N}\left(\mathbf{0},\mathbf{I}{\sigma }_e^2\right)$, where $\mathbf{I}$ is an identity matrix, ${\sigma }_e^2$ is the environmental residual variance.

Simulation model Equation 10 is QTL-based whereas Equation 11 is marker-based prediction model. Simulation and prediction models are also different in how the covariate of additive effects is calculated. The covariate in Equation 10 is independent of the allele frequency in the population while the covariate as elements of ${\mathbf{G}}_{\mathrm{u}}$ in Equation 11 uses the frequency in calculation (Chu et al., 2024). The variances estimated from Equation 11 are statistical parameters whereas the simulated variances in Equation 10 are functional variances (Chu et al., 2024). However, as non-additive genetic effects were not simulated in this example, the functional and statistical variances in example 1 would be identical.

Each founder population (at generation 0) were replicated 5 times, i.e., a total of 50 replicates were simulated. The breeding scheme at generations 1–13 was replicated 10 times per founder population replicate. Variance component estimation was carried out using REML in the DMUAI module of the DMU package (Madsen and Jensen, 2013). The prediction of breeding values was performed with the DMU4 module of the DMU package. Population accuracy of the predicted breeding values were assessed for the individuals in generation 6. The accuracy was the correlation between true $\mathbf{u}$ from the simulated values and predicted $\widehat{\mathbf{u}}$ in (Equation 11). Rate of genetic gain was assessed as the rate of increase in the genetic mean of population from generation 5 to 13, i.e., rate of genetic gain = $\frac{\overline{u_{13}}-\overline{u_5}}{8}$, where $\overline{u_5}$ and $\overline{u_{13}}$ are the genetic means of population at generation 5 and 13.

2.4 Example 2

The investigated factors in example 2 included different prediction models for selection, multi-allelic levels and two levels of ploidy (Table 3). This example used four haplotype founder populations from example 1 for scenarios that had multi-allelic assumption of bi-allelic (100%), and bi-allelic (20%) + quad-allelic (80%) QTL at two ploidy levels of diploids and tetraploids. There were three prediction models investigated, thus in total 2 × 2 × 3 = 12 scenarios simulated. The simulation pipeline of this example is the same as in example 1 and Figure 1. However, the simulation model in example 2 included additive, dominance, and epistatic genetic effects, as in (Equation 9). For simulating epistasis, we assumed $n_{ep}=1000$, with each QTL present in precisely one pair. The functional variance inputs for the simulation model were set for additive ${\sigma }_A^2=1.0$, dominance ${\sigma }_D^2=0.25$, additive × additive ${\sigma }_{AA}^2=0.25$, additive-dominance ${\sigma }_{AD}^2=0.25$, dominance × dominance ${\sigma }_{DD}^2=0.25$, and environmental term ${\sigma }_e^2=2$. The GBLUP models for predicting breeding values included (Equation 11) and two others as follows:

$\mathbf{y}=\mathbf{\text{Xb}}+{\mathbf{Z}}_{\mathbf{u}}\mathbf{u}+{\mathbf{Z}}_{\mathbf{v}}\mathbf{v}+\mathbf{e}(12)$

$\mathbf{y}=\text{Xb}+{\mathbf{Z}}_{\mathbf{u}}\mathbf{u}+{\mathbf{Z}}_{\mathbf{v}}\mathbf{v}+{\mathbf{Z}}_{\text{uv}}\left(\mathbf{\text{uu}}\right)+{\mathbf{Z}}_{\text{uv}}\left(\mathbf{\text{uv}}\right)+{\mathbf{Z}}_{\text{vv}}\left(\mathbf{\text{vv}}\right)+\mathbf{e}(13)$

where $\mathbf{u}$ is the vector of breeding (additive) values as described in Equation 11; $\mathbf{v}$ is the vector of dominance values $\mathbf{v}\sim \mathbf{N}\left(\mathbf{0},{\mathbf{G}}_{\mathrm{v}}{\sigma }_v^2\right)$, where ${\sigma }_v^2$ is dominance variance, and the relationship matrix for dominance (${\mathbf{G}}_{\mathrm{v}}$) is genomic marker-based, which were calculated using AGHmatrix R package (Amadeu et al., 2023). The construction of dominance relationships between individuals used the method by Vitezica et al. (2017) for diploids, and Endelman et al. (2018) for tetraploids. Vector $\left(\mathbf{\text{uu}}\right)$ is additive × additive $\left(\mathbf{\text{uu}}\right)\sim \mathbf{N}\left(\mathbf{0},{\mathbf{G}}_{\text{uu}}{\sigma }_{uu}^2\right)$, where ${\sigma }_{uu}^2$ is additive × additive variance, and the relationship matrix ${\mathbf{G}}_{\text{uu}}={\mathbf{G}}_{\mathrm{u}}\odot {\mathbf{G}}_{\mathrm{u}}$, where ${\mathbf{G}}_{\mathrm{u}}$ is genomic-based relationship as in (Equation 11), $\odot$ is the Hadamard product. Similarly, $\left(\mathbf{\text{uv}}\right)$ is the vector of additive-dominance: $\left(\mathbf{\text{uv}}\right)\sim \mathbf{N}\left(\mathbf{0},{\mathbf{G}}_{\text{uv}}{\sigma }_{uv}^2\right)$, where ${\sigma }_{uv}^2$ is the variance, and ${\mathbf{G}}_{\text{uv}}={\mathbf{G}}_{\mathrm{u}}\odot {\mathbf{G}}_{\mathrm{v}}$. Vector $\left(\mathbf{\text{vv}}\right)$ is dominance × dominance: $\left(\mathbf{\text{vv}}\right)\sim \mathbf{N}\left(\mathbf{0},{\mathbf{G}}_{\text{vv}}{\sigma }_{vv}^2\right)$, where ${\sigma }_{vv}^2$ is the variance, and ${\mathbf{G}}_{\text{vv}}={\mathbf{G}}_{\mathrm{v}}\odot {\mathbf{G}}_{\mathrm{v}}$. Other notations and symbols are the same as in (Equation 11). Variances ${\sigma }_u^2$, ${\sigma }_v^2$, ${\sigma }_{uu}^2$, ${\sigma }_{uv}^2$ and ${\sigma }_{vv}^2$ from Equation 13 are statistical parameters, which are different from functional variances ${\sigma }_A^2$, ${\sigma }_D^2$, ${\sigma }_{AA}^2$, ${\sigma }_{AD}^2$ and ${\sigma }_{DD}^2$ in Equation 9 (Chu et al., 2024).

Similar to example 1, the breeding scheme at generations 1–13 was replicated for 10 times for each of the five founder population replicates. Variance component estimation and prediction of breeding values were carried out by DMUAI and DMU4 module, respectively, of the DMU package (Madsen and Jensen, 2013). Rate of genetic gain was calculated similar as in example 1, except that the total genetic value $g_i$ was used instead of additive values only.

3 Results

Table 4 shows genotypic values of additive ($a$) using (Equation 1) and dominance ($d$) using (Equation 3) for diploids and tetraploids when assuming $n_B=2$ at a locus with allele $B_1$ and $B_2$. With bi-allelic QTL locus of diploids, the genotypic effects of $B_1{B}_1$, $B_1{B}_2$, and $B_2{B}_2$ were $-\left(a_2-a_1\right)$, 0 and $+\left(a_2-a_1\right)$ for $a$, respectively, and 0, $\left(d_1+d_2\right)$ and 0 for $d$, respectively. Table 5 shows epistatic genotypic values for diploids due to interactions between loci $k$ and $l$ when assuming $n_B=2$ at the two loci. The additive-dominance interaction between the pair of loci $k$ and $l$ is the sum of additive × dominance and dominance × additive interaction effects for the two loci.

Table 4

Table 4. Genotypic values of additive ($\mathit{a}$) and dominance ($\mathit{d}$) at one-locus level when assuming bi-allelic QTL.

Table 5

Table 5. Genotypic values of additive × additive $\left(\mathit{a}\mathit{a}\right)$, additive - dominance $\left(\mathit{a}\mathit{d}\right)$, and dominance × dominance $\left(\mathit{d}\mathit{d}\right)$ in two-loci epistatic interactions when assuming bi-allelic QTL.

Table 6 shows rate of genetic gains and accuracy of predicted breeding values in a simplified breeding scheme for diploids and tetraploids when different levels of multiple allelism were assumed in example 1. Surprisingly, accuracy of predicted breeding values was not statistically different for altered levels of QTL multi-allelic assumptions, which occurred in the breeding scheme for both diploids and tetraploids. On the contrary, the rate of genetic gains had increasing tendency with increasing levels of multiple allelic QTL in the scheme for tetraploids whereas the genetic gain did not show this tendency in the scheme for diploids. The increasing tendency in genetic gains was most likely due to increased additive variances with higher multiple allelism in tetraploids. The estimated variances for different scenarios in example 1 could be found in Supplementary Table 1, and the true variances at different generations are in Supplementary Table 2.

Table 6

Table 6. Genetic gain and accuracy in simulation model with additive effects only of example 1.

Table 7 shows the rate of genetic gains in a breeding scheme where different prediction models were used for selection in example 2. In this example, while the simulation model includes additive, dominance, and epistatic interactions between pairs of loci, different prediction models were used for selection. The prediction models with and without non-additive effects did not lead to statistical differences in rate of genetic gains. In about two thirds of replicates, the model could not estimate epistatic effects (Supplementary Table 3). Different multiple allelism did not lead to a significant change in rate of genetic gains. The variance components estimated from different prediction models in example 2 can be found in Supplementary Table 4, and the true variances of total genetic values at different generations are in Supplementary Table 2.

Table 7

Table 7. Genetic gain in simulation model with additive, dominance and epistatic effects of example 2 when different prediction models were used.

4 Discussion

Our simulation models allow the presence of multi-allelic loci, which is a more realistic assumption for QTL variants. Our simulation model for additive genetic effects is basically the sum of allelic effects. When bi-allelic QTL is assumed, the simulation models applying Equations 1, 3, and 5–7 are identical to the assumption in common genetic textbook, e.g., Falconer and Mackay (1996) and AlphaSimR package (Gaynor et al., 2021). For example, AlphaSimR defines effects of genotype $B_1{B}_1$, $B_1{B}_2$ and $B_2{B}_2$ as $–a$, 0 and $+a$ for additive, respectively. These effects would be corresponding to values in Table 4 if $a$ was defined: $a=\left(a_2-a_1\right)$. Here, $\left(a_2-a_1\right)$ is the substitution effect of allele $B_2$ for $B_1$, which is also the definition of $a$ in AlphaSimR and Falconer and Mackay (1996). When bi-allelic QTL is assumed, dominance effect in different ploidy levels is based on a digenic dominance model, which is consistent to Gaynor et al. (2021). Similarly, the epistatic effects in Gaynor et al. (2021) are a special case of our simulation model with bi-allelic assumption.

Interestingly, different levels of multi-allelic assumptions for QTL did not affect accuracy of predicted breeding values based on bi-allelic markers in example 1. This may be due to high density of markers and well-structure population where each clone had many full and half-sibs. Multiple markers can link to a given QTL. Therefore, effects of all alleles at the QTL with multi-allelic state could be estimated using bi-allelic markers. For example, three different bi-allelic marker loci that were closely linked to a QTL could combined to code for up to eight different alleles of the QTL. This could be the reason that regardless of possible multiple alleles in QTL, bi-allelic markers with reasonably high density could predict breeding values in many genomic selection programs (Chu et al., 2019; Samorè and Fontanesi, 2016; Hayes et al., 2013).

In example 1, the differences in rate of genetic gain between different degree of multi-allelisms is primarily due to genetic variances. Although the base population variances are simulated as the same values between two populations, existence of multi-allelism can have higher potential variance, or lower loss of genetic variance under selection. For example, selection led to a removal of a “bad” allele at a QTL in the population. The genetic variance due to QTL would be zero in the bi-allelic case, but might be not in the multi-allelic population. However, maintenance of multi-allelic state might require a bigger effective population size. Otherwise, the multi-allelic state could be lost due to random sampling. This could be the explanation for a higher genetic variance of multi-allelic population in case of tetraploids, but not in diploids in example 1.

While the simulation model including additive, dominance and epistasis was the same for scenarios in example 2, different prediction models (Equations 11–13) were employed for selection. Definition of accuracy of prediction is unclear in literature when different prediction models were used in this case, and particularly when multi-allelic QTL was assumed. Therefore, the rate of genetic gains was used as the main criteria to compare prediction models. Surprisingly, the use of different prediction models did not lead to significant changes in the rate of genetic gains. In other words, the use of correct prediction model for selection, i.e., prediction model and simulation model were more similar, did not improve genetic gains of the breeding scheme. The higher level of multi-allelic assumptions for QTL tended to reduce the genetic gains in example 2, which might be due to lower accuracy of prediction. However, the reduction was not significant.

Nonetheless, examples in this paper are small-scale studies to test ADAM-Multi for multi-allelic features. Many other factors that may affect genetic gain, accuracy of predictions and genetic variances in multi-allelic populations include LD between markers and QTL alleles, population structure, population size, and prediction model. Like other software (Gaynor et al., 2021; Pook et al., 2020; Younis et al., 2023), ADAM-Multi uses functional effects for simulating genotypic values of individuals. Functional effects are independent of allele frequency, thus convenient for studying the consequence of selection in breeding programs (Chu et al., 2024). However, the functional effects and variance parameters cannot be obtained directly by model estimation using real data. Therefore, it is difficult to ensure user-defined parameters for the simulated populations. Just like other software (Gaynor et al., 2021; Pook et al., 2020; Younis et al., 2023), ADAM-Multi is still missing an important feature for a transformation between functional and statistical parameters. More theories are needed to be developed for this transformation, particularly, in case of multi-allelic QTL. Nonetheless, with a more realistic assumption of QTL, ADAM-Multi opens research possibility to study the use of genotyping technology of bi-allelic markers, or the need of new genotyping technology to improve accuracy of selection. Particularly, this assumption of QTL remains very relevant for genomic prediction studies involving multiple breeds and populations. For example, different functional effects of QTL could be assumed in two populations, e.g., González-Diéguez et al. (2021).

This paper presented single-trait models, but the program, ADAM-Multi, can be used for simulating multiple traits with different levels of correlations. For example, the scaling and rescaling procedures in simulation of multiple traits use matrix multiplication, inversion and Cholesky decomposition instead of number multiplication, devision and square root calculations as indicated in this paper. Another note is that the number of alleles $n_B$ could be defined individually for each of QTL (Supplementary Appendix 1). However, this paper assumes a defined $n_B$ for all QTL even when not all alleles in a QTL are segregating. With this assumption, non-segregating alleles’ effects can affect the mean, but the mean can be altered with adding constant values to $\mu$ just like in Equation 9. On the contrary, non-segregating alleles do not affect the functional or statistical variances as the frequencies of these alleles are fixed. In addition, additive effects are centered in our simulation models, but the dominance effects are not. This assumption of dominance leads to a positive effect of heterozygous genotypes, as recommended in Wellmann and Bennewitz (2011).

5 Conclusion

This paper presented a simulation model capable of simulating genotypic effects generalized for multiple allelic models and different ploidy levels. This model accommodates genotypic effects of additive, dominance, and epistasis. When assuming bi-allelic QTL, the generalized model becomes identical to the model assumption in common simulation programs, and in genetic textbooks. This model is integrated in our software ADAM-Multi.

In a small-scale study, we have shown that with a reasonable density of bi-allelic markers and a well-structured population, genomic models can effectively predict breeding values despite the presence of multi-allelic QTL. It was also shown that the inclusion of non-additive genetic effects in the prediction model for selection did not lead to a significant improvement in the rate of genetic gains of a breeding program.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

TC: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Software, Validation, Visualization, Writing–original draft, Writing–review and editing. JJ: Funding acquisition, Methodology, Project administration, Resources, Supervision, Validation, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. The research as a part of project RESPOT was funded by Ministry of Food, Agriculture and Fisheries of Denmark under the Green Development and Demonstrations Program (grant no. 34009-20-1643).

Acknowledgments

We would like to thank Merethe Bagge and the breeders from Danespo for discussions on their breeding program for potatoes, and the presence of multi-allelism in potato population.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2025.1513615/full#supplementary-material

References

Álvarez-Castro, J. M., and Carlborg, O. R. (2007). A unified model for functional and statistical epistasis and its application in quantitative trait loci analysis. Genetics. 176, 1151–1167. doi:10.1534/genetics.106.067348

PubMed Abstract | CrossRef Full Text | Google Scholar

Álvarez-Castro, J. M., and Crujeiras, R. M. (2019). Orthogonal decomposition of the genetic variance for epistatic traits under linkage disequilibrium-applications to the analysis of bateson-dobzhansky-müller incompatibilities and sign epistasis. Front. Genet. 10, 54. doi:10.3389/fgene.2019.00054

PubMed Abstract | CrossRef Full Text | Google Scholar

Álvarez-Castro, J. M., and Yang, R.-C. (2011). Multiallelic models of genetic effects and variance decomposition in non-equilibrium populations. Genetica 139, 1119–1134. doi:10.1007/s10709-011-9614-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Amadeu, R. R., Garcia, A. A. F., Munoz, P. R., and Ferrão, L. F. V. (2023). AGHmatrix: genetic relationship matrices in R. Bioinformatics. 39, btad445. doi:10.1093/bioinformatics/btad445

PubMed Abstract | CrossRef Full Text | Google Scholar

Bengtsson, C., Thomasen, J. R., Kargo, M., Bouquet, A., and Slagboom, M. (2022). Emphasis on resilience in dairy cattle breeding: possibilities and consequences. J. Dairy Sci. 105, 7588–7599. doi:10.3168/jds.2021-21049

PubMed Abstract | CrossRef Full Text | Google Scholar

Biová, J., Kaňovská, I., Chan, Y. O., Immadi, M. S., Joshi, T., Bilyeu, K., et al. (2024). Natural and artificial selection of multiple alleles revealed through genomic analyses. Front. Genet. 14, 1320652. doi:10.3389/fgene.2023.1320652

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C. J., Garrick, D., Fernando, R., Karaman, E., Stricker, C., Keehan, M., et al. (2022). XSim version 2: simulation of modern breeding programs. G3 (Bethesda). 12, jkac032. doi:10.1093/g3journal/jkac032

PubMed Abstract | CrossRef Full Text | Google Scholar

Christensen, O. F., Madsen, P., Nielsen, B., Ostersen, T., and Su, G. (2012). Single-step methods for genomic evaluation in pigs. Animal Int. J. Animal Biosci. 6, 1565–1571. doi:10.1017/s1751731112000742

PubMed Abstract | CrossRef Full Text | Google Scholar

Chu, T. T., Bastiaansen, J. W. M., Berg, P., Romé, H., Marois, D., Henshall, J., et al. (2019). Use of genomic information to exploit genotype-by-environment interactions for body weight of broiler chicken in bio-secure and production environments. Genet. Sel. Evol. 51, 50. doi:10.1186/s12711-019-0493-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Chu, T. T., Henryon, M., Jensen, J., Ask, B., and Christensen, O. F. (2021). Statistical model and testing designs to increase response to selection with constrained inbreeding in genomic breeding programs for pigs affected by social genetic effects. Genet. Sel. Evol. 53, 1. doi:10.1186/s12711-020-00598-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Chu, T. T., Kristensen, P., and Jensen, J. (2024). Simulation of functional additive and non-additive genetic effects using statistical estimates from quantitative genetic models. Hered. (Edinb). 133, 33–42. doi:10.1038/s41437-024-00690-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Chu, T. T., Sørensen, A. C., Lund, M. S., Meier, K., Nielsen, T., and Su, G. (2020). Phenotypically selective genotyping realizes more genetic gains in a rainbow trout breeding program in the presence of genotype-by-environment interactions. Front. Genet. 866, 866. doi:10.3389/fgene.2020.00866

PubMed Abstract | CrossRef Full Text | Google Scholar

Da, Y. (2015). Multi-allelic haplotype model based on genetic partition for genomic prediction and variance component estimation using SNP markers. BMC Genet. 16, 144. doi:10.1186/s12863-015-0301-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Endelman, J. B., Carley, C. A. S., Bethke, P. C., Coombs, J. J., Clough, M. E., da Silva, W. L., et al. (2018). Genetic variance partitioning and genome-wide prediction with allele dosage information in autotetraploid potato. Genetics. 209, 77–87. doi:10.1534/genetics.118.300685

PubMed Abstract | CrossRef Full Text | Google Scholar

Falconer, D. S., and Mackay, T. F. (1996). Introduction to quantitative genetics. England: Pearson Prentice Hall.

Google Scholar

Faux, A. M., Gorjanc, G., Gaynor, R. C., Battagin, M., Edwards, S. M., Wilson, D. L., et al. (2016). AlphaSim: software for breeding program simulation. Plant Genome 9. doi:10.3835/plantgenome2016.02.0013

PubMed Abstract | CrossRef Full Text | Google Scholar

Fisher, R. A. (1930). The genetical theory of natural selection. Oxford: Clarendon.

CrossRef Full Text | Google Scholar

Gaynor, R. C., Gorjanc, G., and Hickey, J. M. (2021). AlphaSimR: an R package for breeding program simulations. G3 (Bethesda). 11, jkaa017. doi:10.1093/g3journal/jkaa017

PubMed Abstract | CrossRef Full Text | Google Scholar

Gebreyesus, G., Sahana, G., Christian Sørensen, A., Lund, M. S., and Su, G. (2020). Novel approach to incorporate information about recessive lethal genes increases the accuracy of genomic prediction for mortality traits. Hered. (Edinb). 125, 155–166. doi:10.1038/s41437-020-0329-5

PubMed Abstract | CrossRef Full Text | Google Scholar

González-Diéguez, D., Legarra, A., Charcosset, A., Moreau, L., Lehermeier, C., Teyssèdre, S., et al. (2021). Genomic prediction of hybrid crops allows disentangling dominance and epistasis. Genetics 218, iyab026. doi:10.1093/genetics/iyab026

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayes, B. J., Lewin, H. A., and Goddard, M. E. (2013). The future of livestock breeding: genomic selection for efficiency, reduced emissions intensity, and adaptation. Trends Genet. 29, 206–214. doi:10.1016/j.tig.2012.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Henryon, M., Liu, H., Berg, P., Su, G., Nielsen, H. M., Gebregiwergis, G. T., et al. (2019). Pedigree relationships to control inbreeding in optimum-contribution selection realise more genetic gain than genomic relationships. Genet. Sel. Evol. 51, 39. doi:10.1186/s12711-019-0475-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Henryon, M., Ostersen, T., Ask, B., Sørensen, A. C., and Berg, P. (2015). Most of the long-term genetic gain from optimum-contribution selection can be realised with restrictions imposed during optimisation. Genet. Sel. Evol. 47, 21. doi:10.1186/s12711-015-0107-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, Y., Chen, S., Wang, X., Liu, M., Iacono, W. G., Hewitt, J. K., et al. (2020). Association analysis and meta-analysis of multi-allelic variants for large-scale sequence data. Genes 11, 586. doi:10.3390/genes11050586

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., Tessema, B. B., Jensen, J., Cericola, F., Andersen, J. R., and Sørensen, A. C. (2019). ADAM-plant: a software for stochastic simulations of plant breeding from molecular to phenotypic level and from simple selection to complex speed breeding programs. Front. Plant Sci. 9, 1926. doi:10.3389/fpls.2018.01926

PubMed Abstract | CrossRef Full Text | Google Scholar

Madsen, P., and Jensen, J. (2013). DMU: a user’s guide. A package for analysing multivariate mixed models. Available at: https://dmu.ghpc.au.dk/dmu/DMU/(Accessed September 12, 2021).

Google Scholar

Massa, A. N., Manrique-Carpintero, N. C., Coombs, J., Haynes, K. G., Bethke, P. C., Brandt, T. L., et al. (2018). Linkage analysis and QTL mapping in a tetraploid russet mapping population of potato. BMC Genet. 19, 87. doi:10.1186/s12863-018-0672-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Pedersen, L., Sørensen, A., Henryon, M., Ansari-Mahyari, S., and Berg, P. (2009). ADAM: a computer program to simulate selective breeding schemes for animals. Livest. Sci. 121, 343–344. doi:10.1016/j.livsci.2008.06.028

CrossRef Full Text | Google Scholar

Pérez-Enciso, M., Ramírez-Ayala, L. C., and Zingaretti, L. M. (2020). SeqBreed: a python tool to evaluate genomic prediction in complex scenarios. Genet. Sel. Evol. GSE 52, 7. doi:10.1186/s12711-020-0530-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Pook, T., Schlather, M., and Simianer, H. (2020). MoBPS - modular breeding program simulator. G3 Genes. Genomes. Genetics 10, 1915–1918. doi:10.1534/g3.120.401193

PubMed Abstract | CrossRef Full Text | Google Scholar

Romé, H., Chu, T. T., Marois, D., Huang, C.-H., Madsen, P., and Jensen, J. (2023). Estimation and consequences of direct-maternal genetic and environmental covariances in models for genetic evaluation in broilers. Genet. Sel. Evol. 55, 58. doi:10.1186/s12711-023-00829-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Samorè, A. B., and Fontanesi, L. (2016). Genomic selection in pigs: state of the art and perspectives. Ital. J. Anim. Sci. 15, 211–232. doi:10.1080/1828051X.2016.1172034

CrossRef Full Text | Google Scholar

Sargolzaei, M., and Schenkel, F. S. (2009). QMSim: a large-scale genome simulator for livestock. Bioinformatics. 25, 680–681. doi:10.1093/bioinformatics/btp045

PubMed Abstract | CrossRef Full Text | Google Scholar

Tessema, B. B., Liu, H., Sørensen, A. C., Andersen, J. R., and Jensen, J. (2020). Strategies using genomic selection to increase genetic gain in breeding programs for wheat. Front. Genet. 11, 578123. doi:10.3389/fgene.2020.578123

PubMed Abstract | CrossRef Full Text | Google Scholar

Thérèse Navarro, A., Tumino, G., Voorrips, R. E., Arens, P., Smulders, M. J. M., van de Weg, E., et al. (2022). Multiallelic models for QTL mapping in diverse polyploid populations. BMC Bioinforma. 23, 67. doi:10.1186/s12859-022-04607-z

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Vitezica, Z. G., Legarra, A., Toro, M. A., and Varona, L. (2017). Orthogonal estimates of variances for additive, dominance, and epistatic effects in populations. Dominance, Epistatic Eff. Populations 206, 1297–1307. doi:10.1534/genetics.116.199406

PubMed Abstract | CrossRef Full Text | Google Scholar

Voorrips, R. E., and Maliepaard, C. A. (2012). The simulation of meiosis in diploid and tetraploid organisms using various genetic models. BMC Bioinforma. 13, 248. doi:10.1186/1471-2105-13-248

PubMed Abstract | CrossRef Full Text | Google Scholar

Weber, S. E., Frisch, M., Snowdon, R. J., and Voss-Fels, K. P. (2023). Haplotype blocks for genomic prediction: a comparative evaluation in multiple crop datasets. Front. plant Sci. 14, 1217589. doi:10.3389/fpls.2023.1217589

PubMed Abstract | CrossRef Full Text | Google Scholar

Wellmann, R., and Bennewitz, J. (2011). The contribution of dominance to the understanding of quantitative genetic variation. Genet. Res. 93, 139–154. doi:10.1017/S0016672310000649

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, R.-C., and Álvarez-Castro, J. M. (2008). Functional and statistical genetic effects with multiple alleles.

Google Scholar

Younis, O. G., Turchetta, M., Ariza Suarez, D., Yates, S., Studer, B., Athanasiadis, I. N., et al. (2023). ChromaX: a fast and scalable breeding program simulator. Bioinformatics 39, btad691. doi:10.1093/bioinformatics/btad691

PubMed Abstract | CrossRef Full Text | Google Scholar

Zaalberg, R. M., Villumsen, T. M., Jensen, J., and Chu, T. T. (2022). Effective selection for lower mortality in organic pigs through selection for total number born and number of dead piglets. Animals. 12, 1796. doi:10.3390/ani12141796

PubMed Abstract | CrossRef Full Text | Google Scholar