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 p l o i d y . 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 = ∑ i B n B t i B a a i B = t 1 a a 1 + t 2 a a 2 + … + 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 = t i B − n p l o i d y 2 2 n p l o i d y (2) where t i B is a raw genotype dosage, or number of copies of allele i B at the locus, n p l o i d y is the ploidy level of genome. Key notations are defined in Table 1.

Similarly, the genotypic model for dominance value ( d ) of a QTL at the locus level is: d = ∑ i B n B t i B d d i B = t 1 d d 1 + t 2 d d 2 + … + 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 n p l o i d y − t i B 2 n p l o i d y 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 .

The simulation model for additive × additive genotypic value a a k l of the two-locus epistatic interaction between the pair of loci k and l : a a k l = t 1 a , l t 2 a , l … t n B a , l ⊗ t 1 a , k t 2 a , k … t n B a , k ⋅ a a 1 k l a a 2 k l … a a n B × n B k l (5) where t 1 a , x t 2 a , x … t n B a , x is a vector of additive covariates for locus x ( k or l ) with n B elements; ⊗ denotes the Kronecker product; denotes the symbol for matrix multiplication; a a 1 k l a a 2 k l … a a n B × n B k l is vector of additive × additive genotypic effects that have n B × 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 a d k l of the epistatic interaction is the sum of additive × dominance and dominance × additive interaction between the pair of loci k and l : a d k l = t 1 d , l t 2 d , l … t n B d , l ⊗ t 1 a , k t 2 a , k … t n B a , k ⋅ a d 1 k l a d 2 k l … a d n B × n B k l + t 1 a , l t 2 a , l … t n B a , l ⊗ t 1 d , k t 2 d , k … t n B d , k ⋅ d a 1 k l d a 2 k l … d a n B × n B k l (6) where t 1 d , x t 2 d , x … t n B d , x is a vector of dominance covariates (Equation 4) for locus x with n B elements; a d 1 k l a d 2 k l … a d n B × n B k l  and  d a 1 k l d a 2 k l … d a n B × n B k l are vectors of additive × dominance ( a d x k l ), and dominance × additive ( d a x k l ) genotypic effects that have n B × n B elements. The value a d x k l is different from d a x k l .

The simulation model for dominance × dominance genotypic value d d k l of the epistatic interaction between the pair of loci k and l : d d k l = t 1 d , l t 2 d , l … t n B d , l ⊗ t 1 d , k t 2 d , k … t n B d , k ⋅ d d 1 k l d d 2 k l … d d n B × n B k l (7) where d d 1 k l d d 2 k l … d d n B × n B k l is vector of dominance × dominance genotypic effects that have n B × 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 = tr t i a ⋅ a + tr t i d ⋅ d + tr t i aa ⋅ aa + tr t i dd ⋅ dd + tr t i ad ⋅ ad + tr t i da ⋅ da (8) where a is a n B × n q t l matrix:

a = a 1 j 1 a 1 j 2 … a 1 n q t l a 2 j 1 a 2 j 2 … a 2 n q t l … … … … a n B j 1 a n B j 2 … a n B n q t l , where n q t l is the number of QTL ;   a i B j q t l is the additive genotypic effect of allele i B at QTL  j q t l ;   t i a  is  a   n q t l × n B matrix:

t i a = t 1 , i a , 1 t 2 , i a , 1 … t n B , i a , 1 t 1 , i a , 2 t 2 , i a , 2 … t n B , i a , 2 … … … … t 1 , i a , n q t l t 2 , i a , n q t l … t n B , i a , n q t l , where t i B , i a , j q t l is the additive covariate of allele i B at QTL  j q t l of individual i .   t i B , i a , j q t l can be calculated based on the genotype of individual i at locus j q t l using Equation 2; tr is the trace of a matrix. Similarly, d is a n B × n q t l matrix of dominance genotypic effects; t i d is a n q t l × n B matrix of dominance covariate of individual i . Matrix aa has a dimension of n B × n B rows and n e p columns:

aa = a a 1 , 1 k l a a 2 , 1 k l … a a n e p , 1 k l a a 1 , 2 k l a a 2 , 2 k l … a a n e p , 2 k l … … … … a a 1 , n B × n B k l a a 2 , n B × n B k l … a a n e p , n B × n B k l , where n e p is number of epistatic interactions between pairs of loci; a a j e p , i B k , i B l k l is the epistatic additive × additive effect at the interaction pair j e p between allele i B k of locus k and allele i B l of locus l .   t i aa is a matrix with n e p rows and n B × n B columns.Row j e p th of t i aa is set up as in Equation 5. For example, t j ep , i aa = t j e p , 1 , i a , l t j e p , 2 , i a , l … t j e p , n B , i a , l ⊗ t j e p , 1 , i a , k t j e p , 2 , i a , k … t j e p , n B , i a , k , where t j e p , 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 e p for individual i .

Similarly, dd is a n B × n B × n e p matrix of dominance × dominance genotypic effects; t i dd  is  a   n e p × n B × n B matrix of dominance × dominance covariates of individual i . The genotypic additive-dominance effects consist of two components: additive × dominance tr t i ad ⋅ ad and dominance × additive tr t i da ⋅ da . Similar to dd , matrices ad and da have dimension of n B × n B × n e p . t i ad and t i da are matrices with n e p rows and n B × n B columns. Row j e p th of t i ad is set up as in Equation 6. For example, t j ep , i ad = t j e p , 1 , i d , l t j e p , 2 , i d , l … t j e p , n B , i d , l ⊗ t j e p , 1 , i a , k t j e p , 2 , i a , k … t j e p , n B , i a , k , where t j e p , 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 e p for individual i . Row j e p th of t i da is: t j ep , i da = t j e p , 1 , i a , l t j e p , 2 , i a , l … t j e p , n B , i a , l ⊗ t j e p , 1 , i d , k t j e p , 2 , i d , k … t j e p , n B , i d , k .

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 a ( n B × n q t l matrix) in ADAM-Multi are:

- Sampling: Each element of matrix a s is sampled from a user-defined normal distribution, e.g., mean of zero and additive variance σ A 2 [or N(0, σ A 2 )]. Matrix a s with the same dimension as a contains starting values of additive effects.

The calculated variances of σ A 2 and σ 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 d ( n B × n q t l matrix) are:

- Sampling: A dominance degree δ i B d , j q t l for allele i B of QTL j q t l is sampled from a user-defined normal distribution N μ δ , σ δ 2 , where μ δ is the mean of dominance degree, and σ δ 2 is the variance of dominance degree. If a user does not provide this distribution, ADAM-Multi default values are N 0.19 , 0.097 as in Wellmann and Bennewitz, (2011). Each element of matrix d * (starting (or prior) values of dominance effects) is generated as: d i B j q t l * = δ i B d , j q t l * a i B j q t l .

Steps for sampling additive × additive effects are.

- Sampling: Each element of matrix aa * [starting (or prior) values of additive × additive effects] is sampled from a user-defined normal distribution N(0, σ A A 2 ).

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 ad * and da * [starting (or prior) values of additive × dominance and dominance × additive effects] is sampled from a user-defined normal distribution N(0, σ A D 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 = μ + g i + e i (9) where g i is the genetic values of individual i that is constructed as in (Equation 8); μ 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.

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.

The simulation model for individuals’ phenotype in example 1 was: y i = μ + ∑ i B n B t i B , i a , 1 a i B j 1 + ∑ i B n B t i B , i a , 2 a i B j 2 + … + ∑ i B n B t i B , i a , n q t l a i B n q t l + e i (10)

Simulation of genetic values in (Equation 10) is an extension of (Equation 1) to the sum effects of n q t l , 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, σ e 2 ) with σ 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 σ A 2 = 1.0 .

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: y = Xb + Z u u + e (11) where y is the vector of individual phenotypes; b is a vector of the fixed effects of individuals’ generation; u is the vector of breeding values u ∼ N 0 , G u σ u 2 , where the relationship matrix for additive genetics ( G u ) was constructed based on bi-allelic markers using method (VanRaden, 2008) for different levels of ploidy assumed in the scenario. The computation of G u was carried out using AGHmatrix R package (Amadeu et al., 2023). X and Z u are design matrices relating individuals to fixed effects and additive genetic effects, respectively. Vector e is an environmental residual term: e ∼ N 0 , I σ e 2 , where I is an identity matrix, σ 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 G 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 u from the simulated values and predicted 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 = u 13 ¯ − u 5 ¯ 8 , where u 5 ¯ and u 13 ¯ are the genetic means of population at generation 5 and 13.

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 e p = 1000 , with each QTL present in precisely one pair. The functional variance inputs for the simulation model were set for additive σ A 2 = 1.0 , dominance σ D 2 = 0.25 , additive × additive σ A A 2 = 0.25 , additive-dominance σ A D 2 = 0.25 , dominance × dominance σ D D 2 = 0.25 , and environmental term σ e 2 = 2 . The GBLUP models for predicting breeding values included (Equation 11) and two others as follows: y = Xb + Z u u + Z v v + e (12) y = Xb + Z u u + Z v v + Z uv uu + Z uv uv + Z vv vv + e (13) where u is the vector of breeding (additive) values as described in Equation 11; v is the vector of dominance values v ∼ N 0 , G v σ v 2 , where σ v 2 is dominance variance, and the relationship matrix for dominance ( G 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 uu is additive × additive uu ∼ N 0 , G uu σ u u 2 , where σ u u 2 is additive × additive variance, and the relationship matrix G uu = G u ⊙ G u , where G u is genomic-based relationship as in (Equation 11), ⊙ is the Hadamard product. Similarly, uv is the vector of additive-dominance: uv ∼ N 0 , G uv σ u v 2 , where σ u v 2 is the variance, and G uv = G u ⊙ G v . Vector vv is dominance × dominance: vv ∼ N 0 , G vv σ v v 2 , where σ v v 2 is the variance, and G vv = G v ⊙ G v . Other notations and symbols are the same as in (Equation 11). Variances σ u 2 , σ v 2 , σ u u 2 , σ u v 2 and σ v v 2 from Equation 13 are statistical parameters, which are different from functional variances σ A 2 , σ D 2 , σ A A 2 , σ A D 2 and σ D D 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.