# The long-term effects of genomic selection: 1. Response to selection, additive genetic variance, and genetic architecture

Yvonne Wientjes, Piter Bijma, Mario Calus, Bas Zwaan, Zulma G. Vitezica,
Joost van den Heuvel

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

## To cite this version:

Yvonne Wientjes, Piter Bijma, Mario Calus, Bas Zwaan, Zulma G. Vitezica, et al.. The long-term effects of genomic selection: 1. Response to selection, additive genetic variance, and genetic architecture. Genetics Selection Evolution, 2022, 54 (1), pp.21. 10.1186/s12711-022-00709-7. hal-03717153

HAL Id: hal-03717153

https://hal.inrae.fr/hal-03717153

Submitted on 8 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# The long-term effects of genomic selection: 1. Response to selection, additive genetic variance, and genetic architecture

Yvonne C. J. Wientjes<fnref n="1" />* ID, Piter Bijma<fnref n="1" />, Mario P. L. Calus<fnref n="1" />, Bas J. Zwaan², Zulma G. Vitezica³ and Joost van den Heuvel²

**Abstract**

**Background:** Genomic selection has revolutionized genetic improvement in animals and plants, but little is known about its long-term effects. Here, we investigated the long-term effects of genomic selection on response to selection, genetic variance, and the genetic architecture of traits using stochastic simulations. We defined the genetic architecture as the set of causal loci underlying each trait, their allele frequencies, and their statistical additive effects. We simulated a livestock population under 50 generations of phenotypic, pedigree, or genomic selection for a single trait, controlled by either only additive, additive and dominance, or additive, dominance, and epistatic effects. The simulated epistasis was based on yeast data.

**Results:** Short-term response was always greatest with genomic selection, while response after 50 generations was greater with phenotypic selection than with genomic selection when epistasis was present, and was always greater than with pedigree selection. This was mainly because loss of genetic variance and of segregating loci was much greater with genomic and pedigree selection than with phenotypic selection. Compared to pedigree selection, selection response was always greater with genomic selection. Pedigree and genomic selection lost a similar amount of genetic variance after 50 generations of selection, but genomic selection maintained more segregating loci, which on average had lower minor allele frequencies than with pedigree selection. Based on this result, genomic selection is expected to better maintain genetic gain after 50 generations than pedigree selection. The amount of change in the genetic architecture of traits was considerable across generations and was similar for genomic and pedigree selection, but slightly less for phenotypic selection. Presence of epistasis resulted in smaller changes in allele frequencies and less fixation of causal loci, but resulted in substantial changes in statistical additive effects across generations.

**Conclusions:** Our results show that genomic selection outperforms pedigree selection in terms of long-term genetic gain, but results in a similar reduction of genetic variance. The genetic architecture of traits changed considerably across generations, especially under selection and when non-additive effects were present. In conclusion, non-additive effects had a substantial impact on the accuracy of selection and long-term response to selection, especially when selection was accurate.

**Background**

Animal breeding has substantially increased the performance of livestock populations over the last century [1, 2]. This has been achieved by selecting the genetically best performing individuals to produce the next generation based on own performance and/or performances of

<footnote>*Correspondence: yvonne.wientjes@wur.nl</footnote>

<footnote n="1">¹ Animal Breeding and Genomics, Wageningen University & Research, 6700 AH Wageningen, The Netherlands
Full list of author information is available at the end of the article</footnote>
relatives. In spite of the strong selection, these pedigree-
based selection methods have proven to be sustainable,
as genetic variation and rates of genetic gain have been
stable for many generations in several animal and plant
species, both in commercial breeding programs and in
experimental selection lines [3–6].

Recently, genomic selection has revolutionized animal
breeding [7, 8]. Within genomic selection, genotypes on
several thousands of DNA markers covering the genome,
along with recorded phenotypes, are used to identify the
genetically best animals. In some breeding programs,
genomic selection has doubled the annual rate of genetic
gain compared to pedigree-based selection [9, 10]. Argu-
ably, genomic selection enables selection for lowly-her-
itable traits [11, 12] and for traits that are difficult or
expensive to measure [13–15], for which pedigree-based
selection is generally not very effective. These properties
have resulted in the rapid implementation of genomic
selection in animal breeding programs worldwide [8,
16–18].

The accuracy of selection and, thereby, the genetic gain
from genomic selection are affected by the genetic archi-
tecture of traits [19–21], i.e. the set of causal loci that
underlie each trait, their allele frequencies, and their sta-
tistical additive effects. The genetic architecture is largely
unknown for most traits, including those under selection
in breeding programmes, but is known to evolve over
time as a result of new mutations and changes in allele
frequencies due to selection and drift [1, 2, 22–26]. When
interactions are present within (dominance) or between
(epistasis) loci, the statistical additive effects (also known
as allele substitution effects) depend on the allele fre-
quencies at the locus itself, as well as on those of inter-
acting loci. This means that functional dominance and
epistatic effects contribute to additive genetic variation,
depending on allele frequencies [27–29]. Although inter-
actions between loci are known and common [30–33],
not much is known about their interaction network or
how those interactions contribute to genetic variance
components or how those contributions change over
generations as a result of drift or selection. To date, the
network of genetic interactions has been most intensively
studied in yeast, where 90% of the loci associated with
a trait were found to be involved in at least one interac-
tion, with only few interactions for most loci, and many
interactions for only a few loci [34–36]. Boone et al. [36]
and Mackay [37] argue that it is likely that this network of
genetic interactions is similar in other species, including
livestock and humans.

We hypothesize that genomic selection accelerates
changes in the genetic architecture of traits across gen-
erations, which can affect long-term genetic gain. The
reason for this is not only that genomic selection is

more effective, but also because the distribution of the
selection pressure across the genome is different with
genomic selection. Classical selection methods based on
pedigree relationships distribute selection pressure rela-
tively evenly across the genome [38]. This is in contrast
to genomic selection, which puts less weight on loci with
rare alleles [38, 39]. Thus, genomic selection methods
select more strongly genomic regions that surround loci
with a large contribution to the additive genetic variance
and may significantly increase changes in allele frequen-
cias at those loci [40]. Therefore, although genomic selec-
tion may substantially accelerate the rate of genetic gain
in the short-term, we expect that by ignoring regions
with a smaller contribution to additive genetic vari-
ance, genomic selection increases the risk of losing rare
favourable alleles or may fail to increase the frequency
of such alleles [41–43]. Loss of rare favourable alleles
reduces genetic variation and genetic gain in the long
term [38] and also limits the potential for future selection
on other traits. However, currently, these expectations
have not been investigated in detail or tested in breeding
populations.

Therefore, the aim of this study was to investigate the
long-term effects of genomic selection on the genetic
architecture of traits. Using simulation, we compared
genomic selection to phenotypic and pedigree-based
selection. We investigated the impact of those selection
methods on the rate of genetic gain, the loss in genetic
variance, and the change in genetic architecture over 50
generations of selection. The results provide more insight
on the long-term evolution of the genetic architecture
and genetic variation of traits under different selection
methods.

Methods

Simulated population

We simulated a livestock population over 50 genera-
tions of selection. As a first step, we constructed a his-
torical population in which selection was absent and
mating was at random, using the QMSim software
[44]. The first 2000 generations (generation – 3050
to – 1050) consisted of 1500 individuals, after which
the size of the population gradually decreased to 100
over 500 generations (generation – 1050 to – 550) to
generate linkage disequilibrium. This was followed by
a gradual increase in population size to 1500 over 500
generations (generation – 550 to – 50). From the last
historical generation (generation – 50), 100 females
and 100 males were randomly sampled and used as
input for our own custom Fortran program, in which
they were randomly mated (mating ratio 1:1) with
a litter size of 10 (5 females and 5 males). In each of
the next 50 discrete generations, 100 females and 100
males were randomly sampled and mated to build-
up mutation-drift equilibrium (generation − 50 to
0), using a randomly selected proportion of 0.2. Gen-
eration 0 formed the base population for the 50 gen-
erations of selection. In the following generations, we
used truncation selection to select the best 100 females
and 100 males, which were randomly mated using a
mating ratio of 1:1 and a litter size of 10 (5 females and
5 males), resulting in a selected proportion of 0.2 for
both females and males. Five selection methods were
used, as explained below. This process was replicated
20 times.

Genome

The simulated genome contained 10 chromosomes of
100 cM each. The number of recombination events per
chromosome was sampled from a Poisson distribution
with on average one recombination per chromosome,
uniformly distributed across the chromosome.

In the historical population, 200,000 randomly
spaced bi-allelic loci per chromosome were simulated
with a recurrent mutation rate of 5 × 10⁻⁵ to maintain
at most two alleles at a locus. The population structure
and mutation rate resulted in a U-shaped allele fre-
quency distribution of the loci in the historical popula-
tion. In the last historical generation, 2000 segregating
loci were randomly selected to become causal loci. By
randomly selecting 200 loci from each of 100 equally-
sized bins based on allele frequency, another set of
20,000 segregating loci were selected as genetic mark-
ers. This resulted in a uniform distribution of the allele
frequency of the markers, reflecting the ascertainment
bias of markers that are typically placed on commercial
marker chips [45–47].

After the historical population, the number of muta-
tions per individual was sampled from a Poisson dis-
tribution with an average of 0.6, which resulted in a
mutational variance of ~0.001σ<sup>2</sup><sub>e</sub> under our simulated
additive model (as explained later), as is often observed
in real populations [48–50]. A random 4000 loci that
did not segregate in the last generation of the histori-
cal population were chosen to be subject to mutation.
The loci and effects of the mutations were recycled to
limit the computational requirement. In each genera-
tion, a locus was drawn from the potential loci that did
not segregate at that time, while maximizing the time
between two mutations at the same locus. As such, each
of the 4000 loci was used on average once each 6 to 7
generations. We believe that recycling the same muta-
tions does not impact the results of our study, because
the vast majority of the mutations are lost in the first
generation due to drift, regardless of their effect.

Genetic and phenotypic values

Three genetic models were used to simulate phenotypic
values; a model with only additive effects (A), a model
with additive and dominance effects (AD), and a model
with additive, dominance, and epistatic effects (ADE).
In the last historical generation, functional (or biolog-
cal) additive and dominance effects were assigned to all
2000 causal loci and to the 4000 loci for mutations. At the
same time, epistatic effects were assigned to 90% of those
loci, as was observed for the yeast data [34].

Functional additive effects (a) were sampled from a
normal distribution with mean 0 and standard devia-
tion 1. Functional dominance effects (d) were simulated
proportional to the additive effect by first sampling a
dominance degree (dd) for each locus from a normal
distribution with mean 0.2 and standard deviation 0.3
[51–53], and then computing the dominance effect of
locus i as $d_i = dd_i|a_i|$. This resulted in mostly positive
dominance effects, with a bit of overdominance, as was
empirically observed in pigs [51].

Only pairwise epistatic effects were simulated, because
higher-order interactions have little effect on the phe-
notypic variance when the allele frequency distribution
is U-shaped [27, 28, 54]. The number of interactions per
locus was sampled using the interaction network found
between the ~6000 genes in yeast [34, 55], with many loci
with few interactions and few loci with many interactions
(Fig. 1). This was done by creating an interaction matrix
from the network in yeast, with elements of 1 when loci
interacted and 0 otherwise. From this matrix, columns
and corresponding rows were selected for all loci with
an interaction. For the interaction between loci *B* and *C*,
nine epistatic degrees (ε) were independently sampled
from a normal distribution with mean 0 and standard
deviation 0.45, one for each of the nine possible two-
locus genotype combinations. The sampled ε were used
to create nine epistatic effects (e) for each interaction
as $e = \epsilon \sqrt{|a_B a_C|}$ (Table 1), resulting in larger epistatic
effects for loci with a larger additive effect. This resulted
in the creation of all types of epistasis, i.e. additive-by-
additive, additive-by-dominance, and dominance-by-
dominance. However, by simulating the epistatic effects
in this random manner, the simulated epistatic effects
also contributed to functional additive or dominance
effects (Table 1). When computing functional additive,
dominance, and epistatic variance components, we first
redistributed the simulated epistatic effects in the correct
underlying functional effects. This was achieved by solv-
ing for each interaction the nine equations in Table 1 for
the eight separate functional additive ($a_B$ and $a_C$), domi-
nance ($d_B$ and $d_C$), additive-by-additive ($k$), additive-by-
dominance ($l$ and $m$) and dominance-by-dominance ($n$)
epistatic effects that were underlying that interaction and
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Histogram of the number of interactions per causal locus</figcaption></figure>

<table><caption>**Table 1** Simulated epistatic model for two-locus interactions</caption><thead><tr><th rowspan="2">Genotype locus <i>B</i></th><th colspan="3">Genotype locus <i>C</i></th></tr><tr><th>CC</th><th>Cc</th><th>cc</th></tr></thead><tbody><tr><td><i>BB</i></td><td><i>e</i><sub>00</sub> = &mu; + <i>a</i><sub>B</sub> + <i>a</i><sub>C</sub> + <i>k</i></td><td><i>e</i><sub>01</sub> = &mu; + <i>a</i><sub>B</sub> + <i>d</i><sub>C</sub> + <i>m</i></td><td><i>e</i><sub>02</sub> = &mu; + <i>a</i><sub>B</sub> &minus; <i>a</i><sub>C</sub> &minus; <i>k</i></td></tr><tr><td><i>Bb</i></td><td><i>e</i><sub>10</sub> = &mu; + <i>d</i><sub>B</sub> + <i>a</i><sub>C</sub> + <i>l</i></td><td><i>e</i><sub>11</sub> = &mu; + <i>d</i><sub>B</sub> + <i>d</i><sub>C</sub> + <i>n</i></td><td><i>e</i><sub>12</sub> = &mu; + <i>d</i><sub>B</sub> &minus; <i>a</i><sub>C</sub> &minus; <i>l</i></td></tr><tr><td><i>bb</i></td><td><i>e</i><sub>20</sub> = &mu; &minus; <i>a</i><sub>B</sub> + <i>a</i><sub>C</sub> &minus; <i>k</i></td><td><i>e</i><sub>21</sub> = &mu; &minus; <i>a</i><sub>B</sub> + <i>d</i><sub>C</sub> &minus; <i>m</i></td><td><i>e</i><sub>22</sub> = &mu; &minus; <i>a</i><sub>B</sub> &minus; <i>a</i><sub>C</sub> + <i>k</i></td></tr></tbody></table>

First, nine epistatic effects ($e_{00}$ to $e_{22}$) were simulated randomly, by sampling for each effect an epistatic degree ($\epsilon$) from a normal distribution and scaling them by the additive effects of the two loci (i.e. $e_{00} = \epsilon_{00} \sqrt{|a_B a_C|}$). Then, those nine epistatic effects were used to estimate the separate functional additive ($a_B$ and $a_C$), dominance ($d_B$ and $d_C$), additive-by-additive ($k$), additive-by-dominance ($l$ and $m$) and dominance-by-dominance ($n$) epistatic effects that were underlying those epistatic effects

adding the additive and dominance effects to the functional additive and dominance effects of the corresponding loci.

The functional genetic effects were combined with the genotypes of the individuals to calculate total genetic values. For each individual, a residual term was also sampled from a normal distribution with mean zero and standard deviation equal to the square root of 1.5 times the variance of total genetic values in the last historical generation, resulting in a broad sense heritability of 0.4 in that generation.

## Statistical effects

The natural and orthogonal interaction approach (NOIA) [56, 57] was applied in each generation to compute statistical additive and dominance effects based on the functional additive, dominance, and epistatic effects of all causal loci (the 2000 segregating causal loci and the 4000 loci for mutations) and their allele frequencies [52]. For each locus *i*, the part of the dominance effect that is statistically additive was calculated as $(1 - 2p_i)d_i$, where $p_i$ is the frequency of the focal allele (i.e. allele *B* for locus *B* in Table 1). For each interaction between loci *B* (with alleles *b* and *B*) and *C* (with alleles *c* and *C*), part of the functional epistasis is converted into statistical additive and statistical dominance effects that were computed from three components: (1) a vector *y* with functional epistatic

effects, $\mathbf{y}' = [e_{00} \ e_{10} \ e_{20} \ e_{01} \ e_{11} \ e_{21} \ e_{02} \ e_{12} \ e_{22}]$, (2) a 9 × 9 diagonal matrix **D** with the expected frequencies of the two-locus haplotypes, assuming that loci segregate independently, and (3) a 9 × 9 matrix **W** with the mean and orthogonal contrasts for the two loci, constructed as $\mathbf{W} = \mathbf{W}_B \otimes \mathbf{W}_C$, with:

$$ \mathbf{W}_B = \begin{bmatrix} 1 & p_{Bb} + 2p_{bb} & \frac{-2p_{Bb}p_{bb}}{p_{BB} + p_{bb} - (p_{BB} - p_{bb})^2} \\ 1 & p_{Bb} + 2p_{bb} - 1 & \frac{4p_{BB}p_{bb}}{p_{BB} + p_{bb} - (p_{BB} - p_{bb})^2} \\ 1 & p_{Bb} + 2p_{bb} - 2 & \frac{-2p_{BB}p_{Bb}}{p_{BB} + p_{bb} - (p_{BB} - p_{bb})^2} \end{bmatrix}, $$

where $p_{BB}$, $p_{Bb}$, and $p_{bb}$ represent the frequencies of the genotypes *BB*, *Bb*, and *bb* for locus *B*. The statistical effects related to the interaction between loci *B* and *C* then follow from:

$$ \mathbf{b}_{BC} = [\mu \ \alpha_{BC}^B \ \delta_{BC}^B \ \alpha_{BC}^C \ (\alpha\alpha)_{BC} \ (\delta\alpha)_{BC} \ \delta_{BC}^C \ (\alpha\delta)_{BC} \ (\delta\delta)_{BC}]' \\ = (\mathbf{W}'\mathbf{D}\mathbf{W})^{-1}\mathbf{W}'\mathbf{D}\mathbf{y}, $$

where $\mu$ is a general mean, $\alpha^x$ is the statistical additive effect related to locus *x*, $\delta^x$ is the statistical dominance effect related to locus *x*, and $\alpha\alpha$, $\alpha\delta$, $\delta\alpha$, and $\delta\delta$ are, respectively, the additive-by-additive, additive-by-dominance, dominance-by-additive, and dominance-by-dominance epistatic effects. Note that the NOIA model was run separately for each pair of interacting loci, such that only the
functional interaction effects were considered and not the functional additive and dominance effects. Therefore $\alpha_{BC}^B = (p_C - q_C)k + 2p_C q_C m + (1 - 2p_B)(p_C - q_C)l + 2p_C q_C (1 - 2p_B)n$, and $\delta_{BC}^B = -(1 - 2p_C)l + 2p_C (1 - p_C)n$; where k, l, m, and n are the additive-by-additive, dominance-by-additive, additive-by-dominance and dominance-by-dominance functional epistatic effects (Table 1), respectively.

The total statistical additive effect at locus i was calculated as:

$$\alpha_i = a_i + (1 - 2p_i)d_i + \sum \alpha_{ij}^i,$$

and the total statistical dominance effect as:

$$\delta_i = d_i + \sum \delta_{ij}^i,$$

where the summations were taken across all interactions that involved locus i.

The statistical additive effect was used to compute the total additive genetic value (i.e. true breeding value) across all loci i of each individual as $A = \sum w_{a_i} \alpha_{i}$, with:

$$w_{a_i} = \begin{cases} p_{Bb} + 2p_{bb} & \text{for genotypes } \{BB, Bb\}. \\ p_{Bb} + 2p_{bb} - 1 & \text{for genotypes } \{Bb, bb\}. \\ p_{Bb} + 2p_{bb} - 2 & \text{for genotypes } \{bb\}. \end{cases}$$

In the same way, the statistical dominance effect was used to compute the total dominance deviation across all loci i of each individual as $D = \sum w_{d_i} \delta_{i}$, with:

$$w_{d_i} = \begin{cases} \frac{-2p_{Bb}p_{bb}}{p_{BB} + p_{bb} - (p_{BB} - p_{bb})^2} & \text{for genotypes } \{BB, Bb\}. \\ \frac{-2p_{BB}p_{Bb}}{p_{BB} + p_{bb} - (p_{BB} - p_{bb})^2} & \text{for genotypes } \{Bb, bb\}. \end{cases}$$

By definition, the variance in A across all individuals is the additive genetic variance, the variance in D across all individuals is the dominance genetic variance, and the variance in total genetic values across all individuals is the total genetic variance. The total genetic variance minus the additive and dominance variance is the epistatic variance.

## Selection methods

Five methods were used to select the sires and dams of the next generation. As a base line for comparison, the first method randomly selected the parents (RAN-DOM) and was meant to capture the impact of drift alone. The second method selected the individuals with the highest phenotypic values to become the parents of the next generation (MASS). The third method selected individuals with the highest estimated breeding values using a pedigree best linear unbiased prediction (BLUP) model that included own performance

information of the selection candidates (PBLUP_OP). The fourth and fifth methods selected individuals with the highest genomic estimated breeding values from a genomic BLUP model that either included own performance information of the selection candidates (GBLUP_OP) or not (GBLUP_NoOP).

Breeding value estimation for the last three methods was performed using the MTG2 software [58]. Each generation, breeding values were estimated simultaneously with estimating the variance components, using the phenotypic information of the previous three generations, and for PBLUP_OP and GBLUP_OP using also phenotypic information of the present generation. The PBLUP method used a relationship matrix based on a pedigree that included all individuals from the present generation and the previous eight generations. The GBLUP methods used a relationship matrix based on marker genotypes of the present generation and the previous three generations, computed using Method 1 of VanRaden [59], with allele frequencies estimated based on the genotype data of those generations. The model for breeding value estimation included a fixed mean, a random additive genetic effect, a random litter effect, and a residual. The random litter effect was included to capture resemblance between full sibs due to non-additive genetic effects, which could otherwise bias the estimated breeding values. Although dominance and epistatic effects were simulated, these were not included in the breeding value estimation model, because additive models are generally used in breeding programs and only the breeding value is transmitted to the offspring.

## Comparing genetic models and selection methods

The three genetic models (A, AD, and ADE) and five selection methods resulted in 15 scenarios that were applied to each of the 20 replicates of the simulated population. The scenarios were compared based on accuracies of selection, phenotypic trend, additive genetic variance, additive genic variance (calculated as the sum of $2p_i(1 - p_i)\alpha_i^2$ across all causal loci i), expected heterozygosity, average minor allele frequency (MAF), and number of segregating causal loci across the 50 generations of selection. Accuracy of selection within a generation was calculated as the correlation between the true and estimated breeding values among animals in that generation.

One of our main aims was to evaluate how fast the genetic architecture of the trait changed due to selection. The genetic architecture can change because: (1) the subset of loci affecting the trait changes due to new mutations and loci becoming fixed, (2) the allele frequencies of those loci change, which can result in changes in the proportion of the additive genetic variance explained by each locus,
or (3) the statistical additive effects of the loci change as a result of changes in allele frequencies and of non-additive effects, which can also change the proportion of additive genetic variance explained by a locus. To quantify changes in genetic architecture, we defined three criteria that each reflected one of those mechanisms, namely: (1) the Jaccard index [60] for the segregating causal loci, (2) the correlation of allele frequencies at those loci between generations, and (3) the correlation of statistical additive effects at those loci between generations. For the first criterion, we calculated the Jaccard index [60] between generation 0 (before selection) and each of the generations after selection as the number of overlapping segregating loci divided by the total number of segregating loci in the two generations. For the second criterion, we calculated the correlation of allele frequencies in generation 0 with those in each subsequent generation, using only the loci that segregated in generation 0 and still segregated in the generation in question. For the third criterion, we calculated the correlation of statistical additive effects in generation 0 with those in a subsequent generation, again including only loci that segregated in both generations.

# Results

## Properties of the simulated population

The distribution of allele frequencies of the segregating causal loci was strongly U-shaped (see Additional file 2: Fig. S1) and comparable to the distribution of allele frequencies of segregating loci that is observed in sequence data of livestock populations [61–64]. In the RANDOM scenario, where no selection was performed, the pattern of allele frequencies remained similar across generations, indicating that the population was approximately in mutation-drift equilibrium. Moreover, the pattern of linkage disequilibrium in the population (see Additional file 2: Fig. S2) was similar to that found in pig and chicken populations [65–67]. This indicates that the effective population size of the simulated population was comparable to that in real livestock populations, which ranges from 40 to 130 [40, 68–70].

With model ADE, epistasis at the functional level was abundant and 49% of the variation in the total genetic value was generated by functional epistatic effects and only 19% by functional additive effects. However, most of the genetic variance at the statistical level was additive (62%) or due to dominance (33%), and only 5% was epistatic variance in generation 0, which is reasonably close to results for litter size in pigs [71]. The broad-sense heritability was set to 0.4 for all genetic models, resulting in a narrow-sense heritability of ~0.25 for model ADE. This heritability was considerably lower than the narrow-sense heritability of ~0.40 for model A and ~0.38 for model AD. Altogether, those parameters indicate that the

genetic architecture that was simulated based on model ADE could represent the genetic architecture of a quantitative trait in a livestock population.

## Accuracy of selection

In the first generation of selection, the accuracy of selection was always highest with genomic selection including own performance (GBLUP_OP) (Fig. 2). The accuracy was ~0.83 for models A and AD, and ~0.72 for model ADE. The lower accuracy for model ADE is a result of the lower narrow-sense heritability for this model. For all genetic models, the accuracy of the pedigree selection scenario with own performance (PBLUP_OP) in generation 1 was ~0.09 lower than with GBLUP_OP, the accuracy of genomic selection without own performance (GBLUP_NoOP) was ~0.13 lower than with GBLUP_OP, and the accuracy of MASS was ~0.21 lower than with GBLUP_OP. As expected, the accuracy of MASS was equal to the square root of the narrow-sense heritability.

Across generations, the accuracy of selection decreased for all scenarios. The decrease was largest in the first generations as a result of the Bulmer effect [72]. Thereafter, the decrease was slightly larger for the genomic selection scenarios (GBLUP_OP and GBLUP_NoOP) than for PBLUP_OP and MASS. As a result, differences in accuracy between the scenarios were smaller after 50 generations of selection than in the first generation. The accuracy decreased fastest under the ADE genetic model, especially for the genomic selection scenarios. Under this genetic model, the accuracies of PBLUP_OP, MASS and GBLUP_OP were similar after 50 generations of selection.

## Genetic gain

Across generations, the average phenotypic value in the population was constant for the RANDOM scenario and increased with selection (Fig. 3). The rates of genetic gain in the first generations resembled the results for accuracy, with the highest values for GBLUP_OP, followed by PBLUP_OP, GBLUP_NoOP, and finally MASS, and smaller values when non-additive effects were present. The rate of genetic gain decreased over generations, but considerably less for MASS than for the other selection methods. Thus, after 50 generations of selection, cumulative genetic gain was greater for MASS than for PBLUP_OP and GBLUP_NoOP under all genetic models, and MASS also outperformed GBLUP_OP under model ADE.

## Additive genetic and genic variance

The additive genetic and genic variances were approximately constant for the RANDOM scenario and decreased with selection (Fig. 4). As expected, the largest
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /></figure>

drop in additive genetic variance was observed in the first generations of selection as a result of the Bulmer effect, similar to what was observed for the accuracy of selection; by more than 20% in the first three generations of selection. The total drop in genetic variance after 50 generations of selection was more or less similar for

GBLUP_OP and GBLUP_NoOP, for which less than 20% of the initial genetic variance was maintained under genetic models A and AD. Under model ADE, more genetic variance was maintained (~24%) for GBLUP_OP and GBLUP_NoOP after 50 generations of selection. Only slightly more genetic variance (~25%) was
maintained with PBLUP_OP, for which the loss in genetic variance was reasonably similar across the three genetic models. With MASS, the loss in genetic variance was considerably less, with ~40% of the variance maintained after 50 generations of selection.

The additive genic variance is not affected by transient effects such as the Bulmer effect [72]. Therefore, the loss in genic variance was smaller than the loss in genetic variance, especially in the first generations (Fig. 4). Except for this difference in the first generations, the trends in additive genic and genetic variance were very similar.

## Number of segregating causal loci

The number of segregating causal loci decreased for the scenarios with selection (Fig. 5). For PBLUP_OP, the number of loci decreased fastest, with a reduction by almost 50% after 50 generations of selection. For GBLUP_OP and GBLUP_NoOP, the decrease was slightly smaller; 42% for GBLUP_OP and 40% for GBLUP_NoOP. For MASS, the decrease was substantially smaller, at only 20%. The loss in segregating loci was slightly smaller when non-additive effects were present. Interestingly, the number of segregating loci in generation 50 was smaller for PBLUP_OP than for GBLUP_OP and GBLUP_NoOP,

while the additive genic variance was slightly larger for PBLUP_OP.

## Average minor allele frequency at segregating causal loci

The additive genic variance depends on the number of segregating causal loci, as well as their MAF. In the first generations of selection, the average MAF of segregating loci increased, especially for PBLUP_OP (Fig. 6). Thereafter, the average MAF decreased and after 50 generations of selection, it was below its initial value, with the smallest values for the GBLUP scenarios. We found that the average MAF of PBLUP_OP and MASS after 50 generations of selection were slightly above the average MAF before selection, but only under the genetic model ADE. The impact of MASS on the average MAF of segregating loci was very small. The higher average MAF for PBLUP_OP can explain the larger additive genic variance for PBLUP_OP than for GBLUP_OP and GBLUP_NoOP, although PBLUP_OP resulted in a smaller number of segregating loci (Fig. 4 vs. Fig. 5).

## Accumulated heterozygosity

In a random mating population, the accumulated heterozygosity depends on the number of segregating causal loci (Fig. 5), their average MAF (Fig. 6) and on

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>
the variation in MAF among those loci (see Additional file 2: Fig. S3 and Additional file 3). As expected, selection resulted in a decrease in the accumulated heterozygosity (Fig. 7). The reduction in accumulated heterozygosity was similar for GBLUP_OP and GBLUP_NoOP, slightly less for PBLUP_OP, and considerably less for MASS. Moreover, the accumulated heterozygosity decreased more slowly when non-additive effects were present. Thus, the decrease in heterozygosity was smaller for pedigree than for genomic selection, and depended on the genetic model.

## Changes in genetic architecture

Across generations, the subset of causal loci underlying the trait (Fig. 8) and their allele frequencies (Fig. 9) and statistical additive effects (Fig. 10) changed. The change in the subset of loci was measured by the Jaccard index, which was substantial, especially in the first generation. Note that approximately 600 new mutations occurred in each generation, most of which were lost immediately. As a result, two consecutive generations already differed in nearly 1200 causal loci. The subset of loci that affected the trait changed considerably with drift (RANDOM) but the change was amplified by selection. After 50 generations, the average Jaccard index was ~0.27 for

RANDOM, ~0.21 for MASS, and between 0.10 and 0.15 for PBLUP_OP, GBLUP_NoOP, and GBLUP_OP. The Jaccard index was slightly higher when non-additive genetic effects were present. Altogether, these results indicate that the subset of loci that affect the trait constantly changed across generations due to new mutations and drift, and that this change was amplified by selection.

Selection also strongly amplified the change in allele frequencies at causal loci compared to drift (Fig. 9) and (see Additional file 4: Figs. S1 to S15). Due to drift alone, the correlation between allele frequencies of causal loci that segregated in both generation 0 and generation 50 was ~0.93 (RANDOM). The change in allele frequencies as a result of selection was largest under model A, with a correlation between the allele frequencies in generations 0 and 50 of only ~0.10 for GBLUP_OP, GBLUP_NoOP, and PBLUP_OP, and of 0.44 for MASS. These correlations were slightly higher under model AD. When epistatic effects were also present, the change in allele frequencies was much smaller, and the correlation was ~0.28 after 50 generations of GBLUP_OP, GBLUP_NoOP, and PBLUP_OP, and 0.66 for MASS.

As a result of the change in allele frequencies, statistical additive effects of the loci changed when non-additive effects were present (Fig. 10) and (see Additional

<figure><img src="image_11.png" /><figcaption>**Fig. 8** Change in the subset of segregating causal loci across generations for the five selection methods and three genetic models. The change in the subset of segregating causal loci is described by the Jaccard index. The five selection methods were: RANDOM selection, MASS selection, PBLUP selection with own performance (PBLUP_OP), GBLUP selection without own performance (GBLUP_NoOP) or with own performance (GBLUP_OP). The three genetic models were a model with only additive effects (A), with additive and dominance effects (AD), or with additive, dominance and epistatic effects (ADE). Results are shown as averages of 20 replicates and the width of the lines represents the average plus and minus one standard error</figcaption></figure>
<figure><img src="image_12.png" /></figure>

file 5: Figs. S1 to S10). These changes were quite limited when only additive and dominance effects were present, with an average correlation of 0.94 between the statistical additive effects in generations 0 and 50 for all selection methods. When epistatic effects were also present, this correlation was much lower, with an average correlation of 0.95 for RANDOM, 0.65 for MASS, 0.51 for PBLUP_OP, 0.47 for GBLUP_NoOP, and 0.45 for GBLUP_OP. Within 10 generations of selection based on GBLUP_OP, GBLUP_NoOP or PBLUP_OP, these correlations had already dropped to ~0.90.

## Discussion

We investigated the long-term effects of genomic selection on the rate of genetic gain, additive genetic variance, and the genetic architecture of quantitative traits. Results showed that, across 50 generations of genomic selection (GBLUP), the accuracy of selection, the rate of genetic gain, the amounts of additive genetic and genic variation, and the number of segregating causal loci decreased. The same trends were also observed for phenotypic (MASS) and pedigree (PBLUP) selection, but the reductions in these parameters were slightly smaller for PBLUP and considerably smaller for MASS. The main results of our

study are summarized in Table 2, which also refers to the most likely mechanisms that underlie the results, which will be further discussed in the following sections.

## Genetic gain

The cumulative genetic gain after 50 generations of GBLUP_OP selection was 8 to 9% higher than with PBLUP_OP, and 16 to 20% higher than with GBLUP_NoOP, mainly as a result of a higher accuracy. Selection resulted in a decrease in the accuracy over generations, which is in agreement with previous research [41, 73]. The drop in accuracy was largest for the model with epistatic effects (as will be further explained later), especially for GBLUP_OP and GBLUP_NoOP. Therefore, the presence of epistatic effects resulted in a larger decrease in rate of genetic gain over generations.

The drop in accuracy over generations was always smaller for MASS than for the other selection methods. Together with the ability of MASS to maintain more genetic variation, this resulted in the highest cumulative genetic gain after 50 generations for MASS when epistasis was present, and almost the highest cumulative gain when epistasis was absent, which agrees with previous research [74, 75].

Based on the causal loci that were segregating after 50 generations of selection, we estimated the theoretical
<figure><img src="image_13.png" /></figure>

maximum genetic gain that could still be achieved when all these loci would become fixed for the favourable allele, using the statistical additive effects of generation 50 (see Additional file 6: Table S1). This theoretical maximum was highest for RANDOM, followed by MASS, and was on average 7.6% and 6.1% higher for GBLUP with or without own performance than for PBLUP with own performance. This suggests that GBLUP is more sustainable in terms of maintaining future genetic gain than PBLUP.

## Genetic variance

All selection methods resulted in a significant loss in genetic variance (Fig. 4). Part of this loss was transient and a result of the Bulmer effect [72]. The reasonably small difference between the genetic and genic variance indicates that this transient loss of genetic variance was limited (see Additional file 2: Fig. S4). Therefore, most of the loss in genetic variance was permanent and resulted from changes in allele frequencies.

Genic variance is a function of the number of segregating causal loci (n), their average heterozygosity ($\bar{H}_E$), the average square of their statistical additive effects ($\bar{\alpha}^2$), and the covariance between their heterozygosity and $\bar{\alpha}^2$.

($Cov(H_E, \alpha^2)$) (see Additional file 3). Although the total loss in genic variance was comparable for GBLUP and PBLUP, PBLUP lost more segregating causal loci than GBLUP (see Additional file 6: Table S2). In contrast, the loss in average heterozygosity level at causal loci was greater for GBLUP than for PBLUP, likely because of stronger family selection with PBLUP, which agrees with its higher level of pedigree inbreeding than with GBLUP (see Additional file 6: Table S3).

Besides the drop in the number of segregating causal loci and their average heterozygosity, the genic variance slightly decreased over generations with GBLUP and PBLUP, as a result of a decrease in $\bar{\alpha}^2$ over generations (see Additional file 6: Table S2). This drop in $\bar{\alpha}^2$ is likely because loci with a larger statistical additive effect were more likely to become fixed over generations, which was stronger when epistasis was present. The covariance between $H_E$ and $\bar{\alpha}^2$ was in general close to zero and contributed only little to the genic variance.

Compared to GBLUP and PBLUP, the loss in genic variance was much smaller for MASS. This was mostly because MASS maintained much more segregating loci, probably because MASS is better able to exploit and maintain rare favorable alleles than GBLUP and PBLUP [41, 76] and because selection pressure on loci is smaller
<table><caption>**Table 2** Summary of the long-term effects of genomic selection and associated mechanisms</caption><thead><tr><td>Long-term effects of genomic selection on</td><td>Simulation results</td><td>Likely or proven mechanisms</td></tr></thead><tbody><tr><td>Rate of genetic gain (Fig. 3)</td><td>Large drop in rate of genetic gain over generations<br/>Epistasis increased the drop in rate of genetic gain over generations</td><td>Loss in additive genetic variance and reduction in accuracy<br/>Epistasis reduced the informativeness of previous generations for breeding value estimation, and increased the level of inbreeding depression</td></tr><tr><td>Loss in additive genetic variance (Fig. 4)</td><td>Higher rate of genetic gain with genomic selection than with pedigree selection<br/>First generations: large drop in additive genetic variance, smaller drop in additive genetic variance<br/>Later generations: large drop in additive genetic and genetic variance</td><td>Accuracy of breeding value estimation is higher with genomic selection than pedigree selection<br/>Bulmer effect, resulting in transient loss in additive genetic variance due to negative covariances between loci<br/>Reduction in the number of segregating loci, because of fixation of alleles due to selection and losing rare favorable alleles as a result of shrinking estimated effects of rare loci towards zero. Moreover, the average heterozygosity level reduced due to selection<br/>Epistasis resulted in fixing a lower number of loci, because the pressure and direction of selection at a locus can change over generations due to changing statistical additive effects</td></tr><tr><td>Change in genetic architecture of traits (Figs. 8, 9, 10)</td><td>Epistasis reduced the loss in additive genetic variance<br/>Similar loss in additive genetic variance with genomic selection than with pedigree selection<br/>Large change in genetic architecture</td><td>Genomic selection maintained more segregating loci, but each of them with a lower MAF than pedigree selection, probably because pedigree selection results in a stronger family selection<br/>Selection changes the allele frequencies of causal loci, thereby changing the subset of segregating causal loci and their statistical additive effects<br/>When epistasis was present, statistical additive effects of causal loci changed across generations, which lowered the change in allele frequency because the pressure and direction of selection at a locus changed across generations<br/>Genomic selection focused more on a subset of genes that change rapidly, but the average change in allele frequency was similar for genomic and pedigree selection</td></tr><tr><td></td><td>Epistasis reduced the change in subset of loci and allele frequencies, but increased the change in statistical additive effects</td><td></td></tr><tr><td></td><td>Subtle differences in change in genetic architecture between pedigree and genomic selection</td><td></td></tr></tbody></table>
with MASS, which reduces the loss of segregating loci as a result of hitchhiking [77].

The loss in genic variance was slightly smaller when non-additive effects were present. With non-additive effects, the statistical additive effects of loci depend on their allele frequencies [37, 78]. For some loci, the sign of the statistical additive effects even changed over generations when epistasis was present (see Additional file 5: Figs. S1 to S10), which changed the direction of selection on these loci and limited the number of loci that became fixed in the population. This resulted in a larger number of segregating loci (Fig. 5) and a higher level of heterozygosity (Fig. 7) after 50 generations of selection when non-additive effects were present (see Additional file 6: Table S2).

## Genetic architecture

Our initial plan was to quantify the change in genetic architecture by the additive genetic correlation between generations. However, this turned out to be very complex, because this correlation depends on the subset of individuals used. For example, the genetic correlation between generations 1 and 10 depends on whether it is estimated based on individuals from generation 1, from generation 10, or both [52]. Therefore, we decided to focus on three measures that reflect the underlying mechanisms, i.e. the changes in the subset of segregating causal loci, in their allele frequencies, and in their statistical additive effects.

Contrary to our expectations and to earlier results [40, 77], the average change in allele frequencies at segregating loci across generations was similar for GBLUP and PBLUP. The variance of the change in allele frequencies at loci was, however, larger for GBLUP than for PBLUP (see Additional file 6: Table S4). These results confirm that GBLUP focusses more on a subset of the genome that changes rapidly in allele frequencies, while PBLUP spreads the selection pressure more evenly across the genome [38, 39].

We hypothesized that the changes in allele frequencies could result in changes in true breeding values over generations when non-additive effects are present. Therefore, we estimated the correlation between true breeding values of individuals from generation 50 with the true breeding values of the same individuals for performance in generations 49 through 47, where true breeding values for a particular generation were calculated using the statistical additive effects of those generations. The correlation was always higher than 0.99 when only additive or additive and dominance effects were present, but substantially lower than 1 for PBLUP and GBLUP when epistasis was present (Table 3; ~0.95 with generation 49, ~0.91 with generation 48, and

<table><caption>**Table 3** Correlation of true breeding values (TBV) of individuals from generation 50 between performance in generation 50 and each of the three previous generations<sup>a</sup> for three genetic models (A, AD, ADE) and five selection methods<sup>b</sup></caption><thead><tr><th rowspan="2"></th><th colspan="3">Correlation in TBV of generation 50 with</th></tr><tr><th>Generation 49</th><th>Generation 48</th><th>Generation 47</th></tr></thead><tbody><tr><td colspan="4"><strong>Model A</strong></td></tr><tr><td>RANDOM</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td></tr><tr><td>MASS</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td></tr><tr><td>PBLUP</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td></tr><tr><td>GBLUP_NoOP</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td></tr><tr><td>GBLUP_OP</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td></tr><tr><td colspan="4"><strong>Model AD</strong></td></tr><tr><td>RANDOM</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td></tr><tr><td>MASS</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td><td>1.00 (0.000)</td></tr><tr><td>PBLUP</td><td>1.00 (0.000)</td><td>0.99 (0.001)</td><td>0.99 (0.001)</td></tr><tr><td>GBLUP_NoOP</td><td>1.00 (0.000)</td><td>0.99 (0.001)</td><td>0.99 (0.001)</td></tr><tr><td>GBLUP_OP</td><td>1.00 (0.000)</td><td>0.99 (0.000)</td><td>0.99 (0.001)</td></tr><tr><td colspan="4"><strong>Model ADE</strong></td></tr><tr><td>RANDOM</td><td>0.99 (0.000)</td><td>0.99 (0.000)</td><td>0.99 (0.000)</td></tr><tr><td>MASS</td><td>0.98 (0.001)</td><td>0.98 (0.002)</td><td>0.97 (0.003)</td></tr><tr><td>PBLUP</td><td>0.95 (0.003)</td><td>0.92 (0.006)</td><td>0.90 (0.007)</td></tr><tr><td>GBLUP_NoOP</td><td>0.96 (0.004)</td><td>0.91 (0.009)</td><td>0.86 (0.012)</td></tr><tr><td>GBLUP_OP</td><td>0.95 (0.003)</td><td>0.91 (0.006)</td><td>0.87 (0.010)</td></tr></tbody></table>

~0.87 with generation 47), with slightly lower values for GBLUP. This indicates that the correlation of true breeding values between generations decreased rapidly, although the correlation of statistical additive effects was very high between subsequent generations (> 0.99, Fig. 10), which is probably because statistical additive effects changed more rapidly for loci that had a high MAF or a large effect (see Additional file 5: Figs. S1 to S10). This phenomenon drastically decreased the informativeness of previous generations for prediction of breeding values, which reduced the accuracy of selection. Thus, recent generations of reference populations are more useful for genomic prediction, not only because they are more closely related to the selection candidates [68, 79, 80], but also because their genetic architecture is more similar to that of the selection candidates. This might explain why it is sometimes beneficial to remove earlier generations from the reference populations [81, 82]. Moreover, the correlation between statistical additive effects cannot be used to investigate

<footnote><sup>a</sup> True breeding values were estimated for individuals from generation 50, for performance in generation 50 and the three previous generations based on the genotypes of the individuals in generation 50 and the statistical additive effects of each of those generations</footnote>

<footnote><sup>b</sup> Results are shown as averages across the 20 replicates with their corresponding standard errors of the mean between brackets</footnote>
the informativeness of previous generations, because a high correlation between statistical additive effects does not necessarily imply that the correlation between true breeding values is high.

## Non-additive effects

A large part of the functional dominance and epistatic effects at causal loci can be converted into additive genetic variance [26, 29, 30, 83]. Due to this conversion, it is difficult to draw conclusions about the magnitude of the functional epistasis in a population based on the level of epistatic variance [84, 85]. In our study, for example, the epistatic variance only explained 5% of the genetic variance, although almost 50% of the variation in the total genetic value was generated by functional epistatic effects.

The conversion of non-additive effects into statistical additive effects depends on allele frequencies, with a larger proportion of the non-additive effects converted into statistical additive effects when allele frequencies are closer to 0 or 1 [27–29]. As a result, the MAF of a locus was negatively correlated with the absolute value of its statistical additive effect in our simulations, already before selection, although functional effects were simulated independently of allele frequencies (Fig. 11). A negative correlation between MAF and the size of the estimated additive effects of loci is often observed in empirical studies [86–88], and our results show that the presence of non-additive effects can contribute to explaining this finding and that it is not necessarily a result of selection.

Little is known about the structure and network of epistatic interactions. We only simulated pairwise interactions and mimicked the genetic interaction network observed in yeast, with many loci that have few interactions and few loci that have many interactions. Although studied in less detail than in yeast, similar interaction networks have been reported for *C. elegans* [89], *Drosophila* [31] and mice [90], and are also found for proteins [35], which led Boone et al. [36] and Mackay [37] to argue that it is likely that the interaction network between genetic loci is similar in other species such as livestock and humans. Thus, we used the knowledge that is available on interaction effects from yeast [34] to make the simulations as realistic as possible. However, we still had to make simplifications such as only including pair-wise interactions and independence of interaction effects between pairs. We expect that the trend in genetic gain, genetic variance, and genetic architecture for the three genetic models is likely correct because it follows expectations, however, the magnitude of the differences between models and selection methods could be affected by those simplifications.

## Genetic evaluation methods

To estimate genomic relationships, we used the allele frequencies for the individuals that were included in the genomic relationship matrix, i.e. the selection candidates and the previous three generations. In an additive analysis model, the use of different allele frequencies affects the estimates of the variance components but not the ranking of the estimated breeding values [91, 92] and, therefore, not the response to selection. The

<figure><img src="image_14.png" /><figcaption>**Fig. 11** Correlation between the absolute statistical additive effect and the minor allele frequency at causal loci for the three genetic models. The three genetic models were a model with only additive effects (A), with additive and dominance effects (AD), or with additive, dominance and epistatic effects (ADE)</figcaption></figure>
allele frequencies of the markers followed a uniform
distribution, as is common for loci that are included on
most commercial marker chips [45–47], which is differ-
ent from the distribution of the causal loci, which was
U-shaped. This difference could result in a slight bias in
the estimated variance components and breeding values
for the GBLUP scenarios [93, 94], which could have had a
slightly negative effect on the long-term response for the
GBLUP scenarios.

The models used in our study to estimate breeding val-
ues included only an additive effect. The availability of
genomic data enables inclusion of non-additive effects
in the breeding value estimation model [95–97]. We
tested the benefit of including a dominance effect in the
GBLUP model with own performance for the genetic
models AD and ADE. With the AD genetic model, inclu-
sion of a dominance effect in the GBLUP model had a
negligible effect on genetic gain and loss of genetic vari-
ance (see Additional file 2: Figs. S5 and S6). With the
ADE genetic model, inclusion of a dominance effect in
GBLUP resulted in slightly higher long-term response to
selection, which became apparent after about 20 genera-
tions. This is likely related to the much larger dominance
variance in the ADE model compared to the AD model,
which probably makes it more difficult for the random
litter effect to capture all dominance genetic variance,
resulting in a slight bias in estimated breeding values
when a dominance effect is not included, which accumu-
lates over generations. With the ADE genetic model, the
cumulative genetic gain for the GBLUP model with dom-
inance effects was comparable to MASS selection after 50
generations of selection. This result shows that inclusion
of a dominance effect in the breeding value estimation
model can be beneficial in the long-term when non-addi-
tive effects are large.

Relevance for breeding programs

Our results show a larger loss in genetic variance than
generally observed in actual populations, which is often
observed in simulation studies e.g. [41, 43, 74, 76]. The
loss in genetic variance was mainly due to a reduction in
the number of segregating causal loci and in the mean
heterozygosity at these loci. These results suggest that the
average change in allele frequency was larger in our sim-
ulations than in actual populations. Moreover, the large
differences between the RANDOM scenario and the
selected scenarios suggest that the large changes in allele
frequencies were mainly due to selection, rather than
drift. The change in allele frequency at a locus due to
selection depends on the proportion of genetic variance
explained by the locus and on the level of linkage dise-
quilibrium among loci that are under selection. Hence, in
our simulations, a typical causal locus may have explained

more variation than in actual populations, due to either a
larger effect or a higher MAF, and the smaller genome of
roughly 1/3 of the size of a typical livestock genome may
have increased the effect of linkage. At present, we have
insufficient knowledge of the number of loci and the joint
distribution of the allele effects and allele frequencies in
livestock populations to draw strong conclusions. It is
also difficult to predict whether the larger loss in genetic
variance than observed in actual populations affects the
comparison of the selection methods. To be more in line
with reality, a larger genome with more causal loci and
lower MAF could be simulated, resulting in a less severe
reduction in genetic variance. However, an agreement of
the trend in genetic variance between simulations and
actual populations still does not prove that the changes in
genetic architecture observed in simulation match those
of actual populations.

In our simulations, MASS outperformed GBLUP after
50 generations when both dominance and epistatic
effects were present. Note, however, that we compared
selection responses per generation, while the generation
interval may differ between selection methods. Especially
for GBLUP without own performance, the generation
interval can be substantially reduced in some livestock
populations, such as dairy cattle [9, 10]. Moreover, com-
mercial breeding programs typically control the rate of
inbreeding, for example by optimal contribution selec-
tion [98–101], which limits the loss in long-term genetic
gain.

The focus of this study was not on optimization of
breeding programs for long-term genetic gain, which is
not realistic in practice because of competition between
breeding companies, and optimization of long-term gain
requires sacrificing short-term gain. Instead, our focus
was on the long-term consequences of genomic selection
compared to selection strategies that have a longer his-
tory in livestock populations, such as MASS and PBLUP,
which have been proven to be sustainable for a relatively
large number of generations.

Conclusions

Our results show that short-term response was always
greatest with GBLUP, while long-term response was
greater for MASS than for GBLUP when epistasis was
present, and was always greater for MASS than for
PBLUP. This was mainly the result of a much larger loss
in genetic variance and number of segregating loci with
GBLUP and PBLUP than with MASS. The genetic gain
of PBLUP with own performance records was always in
between that of GBLUP with and without own perfor-
mance records. GBLUP and PBLUP showed a similar loss
in genetic variance, but the underlying mechanism was
different with GBLUP maintaining more loci but with a
lower MAF than PBLUP. The maximum genetic gain that could still be obtained after 50 generations was higher for GBLUP selection than for PBLUP, which suggests that GBLUP maintains long-term genetic gain better than PBLUP. Changes in the genetic architecture of the trait, i.e. Jaccard index of segregating causal loci, correlation in allele frequencies, and correlations in statistical additive effects across generations, were strongly amplified by selection but, in contrast to our hypothesis, comparable for GBLUP and PBLUP. Non-additive effects were relatively unimportant in the short-term but had a substantial impact on the accuracy and genetic gain when multiple generations were included in the reference population and selection was accurate.

## Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1186/s12711-022-00709-7.

**Additional file 1.** Programs and seeds to simulate data. This file contains the QMSim input file, Fortran programs and seeds used to select the markers and causal loci, to simulate functional effects and genotypes and phenotypic values of new generations, and the interaction matrix used to simulate epistatic effects.

**Additional file 2: Figure S1.** Allele frequency distribution of segregating causal loci over 50 generations without selection. **Figure S2.** Extent of LD ($r^2$) in the simulated population before selection as a function of distance in one random replicate. **Figure S3.** Trend in the variation in minor allele frequency (MAF) of segregating causal loci for the five selection methods and three genetic models. The five selection methods were: RANDOM selection, MASS selection, PBLUP selection with own performance (PBLUP_OP), GBLUP selection without own performance (GBLUP_NoOP) or with own performance (GBLUP_OP). The three genetic models were a model with only additive effects (A), with additive and dominance effects (AD), or with additive, dominance and epistatic effects (ADE). Results are shown as averages of 20 replicates and the width of the lines represents the average plus and minus one standard error. **Figure S4.** The difference between additive genetic and additive genetic variance which represents a transient loss in genetic variance for the five selection methods and three genetic models. The five selection methods were: RANDOM selection, MASS selection, PBLUP selection with own performance (PBLUP_OP), GBLUP selection without own performance (GBLUP_NoOP) or with own performance (GBLUP_OP). The three genetic models were a model with only additive effects (A), with additive and dominance effects (AD), or with additive, dominance and epistatic effects (ADE). Results are shown as averages of 20 replicates and the width of the lines represents the average plus and minus one standard error. **Figure S5.** Phenotypic trend for the GBLUP model with own performance records and with and without a dominance effect for the genetic models with non-additive effects. The phenotypic trend is scaled by the additive genetic standard deviation in the generation before selection in order to make the results comparable across the genetic models. The two genetic models were a model with additive and dominance effects (AD), or with additive, dominance and epistatic effects (ADE). Results are shown as averages of 20 replicates and the width of the lines represents the average plus and minus one standard error. **Figure S6.** Trend in additive genetic (A, B) and additive genetic (C, D) variance for the GBLUP model with own performance and with and without a dominance effect for the genetic models with non-additive effects. The trend is scaled by the additive genetic or additive genetic variance in the generation before selection in order to make the results comparable across the genetic models. The two genetic models were a model with additive and dominance effects (AD), or with additive, dominance and epistatic effects (ADE). Results are shown as averages of 20 replicates and

the width of the lines represents the average plus and minus one standard error.

**Additional file 3.** Decomposition of additive genetic variance. This file provides a theoretical decomposition of the additive genetic variance.

**Additional file 4: Figure S1.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive effects (Model A) under RANDOM selection. **Figure S2.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive effects (Model A) under MASS selection. **Figure S3.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive effects (Model A) under PBLUP selection with own performance (PBLUP_OP). **Figure S4.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive effects (Model A) under GBLUP selection without own performance (GBLUP_NoOP). **Figure S5.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive and dominance effects (Model AD) under RANDOM selection. **Figure S7.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive and dominance effects (Model AD) under MASS selection. **Figure S8.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive and dominance effects (Model AD) under PBLUP selection with own performance (PBLUP_OP). **Figure S9.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive and dominance effects (Model AD) under GBLUP selection without own performance (GBLUP_OP). **Figure S10.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive and dominance effects (Model AD) under GBLUP selection with own performance (GBLUP_OP). **Figure S11.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive, dominance and epistatic effects (Model ADE) under RANDOM selection. **Figure S12.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive, dominance and epistatic effects (Model ADE) under MASS selection. **Figure S13.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive, dominance and epistatic effects (Model ADE) under PBLUP selection with own performance (PBLUP_OP). **Figure S14.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive, dominance and epistatic effects (Model ADE) under GBLUP selection without own performance (GBLUP_NoOP). **Figure S15.** Scatterplot of allele frequencies of all causal variants in different generations for the genetic model with additive, dominance and epistatic effects (Model ADE) under GBLUP selection with own performance (GBLUP_OP).

**Additional file 5: Figure S1.** Scatterplot of statistical additive effects in different generations for the genetic model with additive and dominance effects (Model AD) under RANDOM selection. **Figure S2.** Scatterplot of statistical additive effects in different generations for the genetic model with additive and dominance effects (Model AD) under MASS selection. **Figure S3.** Scatterplot of statistical additive effects in different generations for the genetic model with additive and dominance effects (Model AD) under PBLUP selection with own performance (PBLUP_OP). **Figure S4.** Scatterplot of statistical additive effects in different generations for the genetic model with additive and dominance effects (Model AD) under GBLUP selection without own performance (GBLUP_NoOP). **Figure S5.** Scatterplot of statistical additive effects in different generations for the genetic model with additive and dominance effects (Model AD) under GBLUP selection with own performance (GBLUP_OP). **Figure S6.** Scatterplot of statistical additive effects in different generations for the genetic model with additive, dominance and epistatic effects (Model ADE) under RANDOM selection. **Figure S7.** Scatterplot of statistical additive effects in different generations for the genetic model with additive, dominance and epistatic effects (Model ADE) under MASS selection. **Figure S8.** Scatterplot of statistical additive effects in different generations for the genetic model with additive, dominance and epistatic effects (Model ADE) under PBLUP selection with own performance (PBLUP_OP).
selection with own performance (PBLUP_OP). **Figure S9.** Scatterplot of statistical additive effects in different generations for the genetic model with additive, dominance and epistatic effects (Model ADE) under GBLUP selection without own performance (GBLUP_NoOP). **Figure S10.** Scatterplot of statistical additive effects in different generations for the genetic model with additive, dominance and epistatic effects (Model ADE) under GBLUP selection with own performance (GBLUP_OP).

**Additional file 6: Table S1.** Maximum genetic gain¹ that is still possible after 50 generations of selection for the five selection methods and three genetic models². The five selection methods were: RANDOM selection, MASS selection, PBLUP selection with own performance (PBLUP_OP), GBLUP selection without own performance (GBLUP_NoOP) or with own performance (GBLUP_OP). The three genetic models were a model with only additive effects (A), with additive and dominance effects (AD), or with additive, dominance and epistatic effects (ADE). ¹The maximum genetic gain in generation 50 is estimated as the genetic gain when all loci would be fixed for the favourable allele, using the statistical additive effects of generation 50 and neglecting mutations. ²Results are shown as averages across the 20 replicates with their corresponding standard errors of the mean between brackets. **Table S2.** Percentual change in the components of the genetic variance after 10 and 50 generations of selection for the five selection methods and three genetic models¹. The five selection methods were: RANDOM selection, MASS selection, PBLUP selection with own performance (PBLUP_OP), GBLUP selection without own performance (GBLUP_NoOP) or with own performance (GBLUP_OP). The three genetic models were a model with only additive effects (A), with additive and dominance effects (AD), or with additive, dominance and epistatic effects (ADE). ¹Results are shown as averages of 20 replicates with their corresponding standard errors of the mean between brackets. Increases in the value of a component are represented in bold. **Table S3.** Average pedigree inbreeding coefficient after 50 generations of selection for the five selection methods and three genetic models¹. The five selection methods were: RANDOM selection, MASS selection, PBLUP selection with own performance (PBLUP_OP), GBLUP selection without own performance (GBLUP_NoOP) or with own performance (GBLUP_OP). The three genetic models were a model with only additive effects (A), with additive and dominance effects (AD), or with additive, dominance and epistatic effects (ADE). ¹Results are shown as averages of 20 replicates with their corresponding standard errors of the mean between brackets.

**Table S4.** Average and variance of change in allele frequency of causal loci¹ across 50 generations of selection for the five selection methods and three genetic models². The five selection methods were: RANDOM selection, MASS selection, PBLUP selection with own performance (PBLUP_OP), GBLUP selection without own performance (GBLUP_NoOP) or with own performance (GBLUP_OP). The three genetic models were a model with only additive effects (A), with additive and dominance effects (AD), or with additive, dominance and epistatic effects (ADE). ¹Causal loci included only the causal loci segregating in generation 0. ²Results are shown as averages of 20 replicates with their corresponding standard errors of the mean between brackets.

## Acknowledgements

Not applicable.

## Authors' contributions

YCWJ obtained funding for this study. YCWJ, PB, MPLC, BJZ, ZGV and JH (all authors) participated in the design of the study. YCWJ performed the simulations and statistical analyses and wrote the first draft of the paper. YCWJ, PB, MPLC, BJZ, ZGV and JH were involved in the interpretation of the results. All authors read and approved the final manuscript.

## Funding

This publication is part of the project 'Revolution of traits? Quantifying the genetic change in traits over generations as a result of Genomic Selection' (with project number 16774) of the research programme Veni which is (partly) financed by the Dutch Research Council (NWO). The use of the HPC cluster has been made possible by CAT-AgroFood (Shared Research Facilities Wageningen UR).

## Availability of data and materials

All scripts used to generate the data during this study are included in Additional file 1. This file contains the QMSim input file, Fortran programs and seeds used to select the markers and causal loci, to simulate functional effects and genotypes and phenotypic values of new generations, and the interaction matrix used to simulate epistatic effects.

## Declarations

### Ethics approval and consent to participate

Not applicable.

### Consent for publication

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

### Author details

¹Animal Breeding and Genomics, Wageningen University & Research, 6700 AH Wageningen, The Netherlands. ²Laboratory of Genetics, Wageningen University & Research, 6700 AH Wageningen, The Netherlands. ³INRAE, INP, UMR 1388 GenPhySE, 31326 Castanet-Tolosan, France.

Received: 9 April 2021 Accepted: 10 February 2022

Published online: 07 March 2022

## References

1. Hill WG. Is continued genetic improvement of livestock sustainable? Genetics. 2016;202:877-81.

2. Hill WG, Kirkpatrick M. What animal breeding has taught us about evolution. Annu Rev Ecol Evol Syst. 2010;41:1-19.

3. Havenstein GB, Ferket PR, Qureshi MA. Growth, livability, and feed conversion of 1957 versus 2001 broilers when fed representative 1957 and 2001 broiler diets. Poult Sci. 2003;82:1500-8.

4. Dudley JW, Lambert RJ. 100 generations of selection for oil and protein in corn. Plant Breed Rev. 2003;24:79-110.

5. Havenstein GB, Ferket PR, Qureshi MA. Carcass composition and yield of 1957 versus 2001 broilers when fed representative 1957 and 2001 broiler diets. Poult Sci. 2003;82:1509-18.

6. Beniwal BK, Hastings IM, Thompson R, Hill WG. Estimation of changes in genetic parameters in selected lines of mice using REML with an animal model. 1. Lean mass. Heredity (Edinb). 1992;69:352-60.

7. Meuwissen THE, Hayes BJ, Goddard ME. Prediction of total genetic value using genome-wide dense marker maps. Genetics. 2001;157:1819-29.

8. Meuwissen THE, Hayes BJ, Goddard ME. Genomic selection: a paradigm shift in animal breeding. Anim Front. 2016;6:6-14.

9. García-Ruiz A, Cole JB, VanRaden PM, Wiggans GR, Ruiz-López FJ, Van Tassell CP. Changes in genetic selection differentials and generation intervals in US Holstein dairy cattle as a result of genomic selection. Proc Natl Acad Sci USA. 2016;113:E3995-4004.

10. Schaeffer L. Strategy for applying genome-wide selection in dairy cattle. J Anim Breed Genet. 2006;123:218-23.

11. Wolc A, Stricker C, Arango J, Settar P, Fulton JE, O'Sullivan NP, et al. Breeding value prediction for production traits in layer chickens using pedigree or genomic relationships in a reduced animal model. Genet Sel Evol. 2011;43:5.

12. Calus MPL, Meuwissen THE, de Roos APW, Veerkamp RF. Accuracy of genomic selection using different methods to define haplotypes. Genetics. 2008;178:553-61.

13. Goddard ME, Hayes BJ. Mapping genes for complex traits in domestic animals and their use in breeding programmes. Nat Rev Genet. 2009;10:381-91.

14. Daetwyler HD, Swan AA, van der Werf JH, Hayes BJ. Accuracy of pedigree and genomic predictions of carcass and novel meat quality traits in multi-breed sheep data assessed by cross-validation. Genet Sel Evol. 2012;44:33.
15. Calus MPL, de Haas Y, Pszczola M, Veerkamp RF. Predicted accuracy of and response to genomic selection for new traits in dairy cattle. *Animal. 2013;7:183–91.

16. Knol EF, Nielsen B, Knap PW. Genomic selection in commercial pig breeding. *Anim Front. 2016;6:15–22.

17. Wolc A, Kranis A, Arango J, Settar P, Fulton JE, O'Sullivan NP, et al. Implementation of genomic selection in the poultry industry. *Anim Front. 2016;6:23–31.

18. Hayes BJ, Bowman PJ, Chamberlain AJ, Goddard ME. Invited review: genomic selection in dairy cattle: progress and challenges. *J Dairy Sci. 2009;92:433–43.

19. Daetwyler HD, Pong-Wong R, Villanueva B, Woolliams JA. The impact of genetic architecture on genome-wide evaluation methods. *Genetics. 2010;185:1021–31.

20. Wientjes YCJ, Calus MPL, Goddard ME, Hayes BJ. Impact of QTL properties on the accuracy of multi-breed genomic prediction. *Genet Sel Evol. 2015;47:42.

21. Hayes BJ, Pryce JE, Chamberlain AJ, Bowman PJ, Goddard ME. Genetic architecture of complex traits and accuracy of genomic prediction: coat colour, milk-fat percentage, and type in Holstein cattle as contrasting model traits. *PLoS Genet. 2010;6:e1001139.

22. Wright S. Evolution in Mendelian populations. *Genetics. 1931;16:97–159.

23. Robertson A. A theory of limits in artificial selection. *Proc R Soc Lond B Biol Sci. 1960;153:234–49.

24. Falconer DS, Mackay TFC. Introduction to quantitative genetics. 4th ed. Harlow: Pearson Education Limited; 1996.

25. Hansen TF, Álvarez-Castro JM, Carter AJR, Hermisson J, Wagner GP. Evolution of genetic architecture under directional selection. *Evolution. 2006;60:1523–36.

26. Le Rouzic A, Carlborg Ö. Evolutionary potential of hidden genetic variation. *Trends Ecol Evol. 2008;23:33–7.

27. Hill WG, Goddard ME, Visscher PM. Data and theory point to mainly additive genetic variance for complex traits. *PLoS Genet. 2008;4:e1000008.

28. Mäki-Tanila A, Hill WG. Influence of gene interaction on complex trait variation with multilocus models. *Genetics. 2014;198:355–67.

29. Barton NH, Turelli M. Effects of genetic drift on variance components under a general model of epistasis. *Evolution. 2004;58:2111–32.

30. Carlborg Ö, Jacobsson L, Ahgren P, Siegel P, Andersson L. Epistasis and the release of genetic variation during long-term selection. *Nat Genet. 2006;38:418–20.

31. Huang W, Richards S, Carbone MA, Zhu D, Anholt RRH, Ayroles JF, et al. Epistasis dominates the genetic architecture of Drosophila quantitative traits. *Proc Nat Acad Sci USA. 2012;109:15553–9.

32. Flint J, Mackay TFC. Genetic architecture of quantitative traits in mice, flies, and humans. *Genome Res. 2009;19:723–33.

33. Carlborg Ö, Haley CS. Epistasis: too often neglected in complex trait studies? *Nat Rev Genet. 2004;5:618–25.

34. Costanzo M, VanderSluis B, Koch EN, Baryshnikova A, Pons C, Tan G, et al. A global genetic interaction network maps a wiring diagram of cellular function. *Science. 2016;353:aaf1420.

35. Tong AHY, Lesage G, Bader GD, Ding H, Xu H, Xin X, et al. Global mapping of the yeast genetic interaction network. *Science. 2004;303:808–13.

36. Boone C, Bussey H, Andrews BJ. Exploring genetic interactions and networks with yeast. *Nat Rev Genet. 2007;8:437–49.

37. Mackay TFC. Epistasis and quantitative traits: using model organisms to study gene–gene interactions. *Nat Rev Genet. 2014;15:22–33.

38. Goddard ME. Genomic selection: prediction of accuracy and maximisation of long term response. *Genetica. 2009;136:245–57.

39. Bijma P. Long-term genomic improvement—new challenges for population genetics. *J Anim Breed Genet. 2012;129:1–2.

40. Heidaritabar M, Vereijken A, Muir WM, Meuwissen T, Cheng H, Megens H-J, et al. Systematic differences in the response of genetic variation to pedigree and genome-based selection methods. *Heredity (Edinb). 2014;113:503–13.

41. Jannink J-L. Dynamics of long-term genomic selection. *Genet Sel Evol. 2010;42:35.

42. De Beukelaer H, Badke Y, Fack V, De Meyer G. Moving beyond managing realized genomic relationship in long-term genomic selection. *Genetics. 2017;206:1127–38.

43. Liu H, Meuwissen THE, Sørensen AC, Berg P. Upweighting rare favourable alleles increases long-term genetic gain in genomic selection programs. *Genet Sel Evol. 2015;47:19.

44. Sargolzaei M, Schenkel FS. QMSim: a large-scale genome simulator for livestock. *Bioinformatics. 2009;25:680–1.

45. Groenen MAM, Megens H-J, Zare Y, Warren WC, Hillier LW, Crooijmans RPMA, et al. The development and characterization of a 60K SNP chip for chicken. *BMC Genomics. 2011;12:274.

46. Matukumalli LK, Lawley CT, Schnabel RD, Taylor JF, Allan MF, Heaton MP, et al. Development and characterization of a high density SNP genotyping assay for cattle. *PLoS One. 2009;4:e5350.

47. Ramos AM, Crooijmans RPMA, Affara NA, Amaral AJ, Archibald AL, Beever JE, et al. Design of a high density SNP genotyping assay in the pig using SNPs identified and characterized by next generation sequencing technology. *PLoS One. 2009;4:e6524.

48. Lynch M, Walsh B. Genetics and analysis of quantitative traits. *Sunderland: Sinauer Associates Inc.; 1998.

49. Hill WG. Predictions of response to artificial selection from new mutations. *Genet Res. 1982;40:255–78.

50. Houle D, Morikawa B, Lynch M. Comparing mutational variabilities. *Genetics. 1996;143:1467–83.

51. Bennewitz J, Meuwissen THE. The distribution of QTL additive and dominance effects in porcine F2 crosses. *J Anim Breed Genet. 2010;127:171–9.

52. Duenk P, Bijma P, Calus MPL, Wientjes YCJ, van der Werf HJH. The impact of non-additive effects on the genetic correlation between populations. *G3 (Bethesda). 2020;10:783–95.

53. Sun X, Mumm RH. Method to represent the distribution of QTL additive and dominance effects associated with quantitative traits in computer simulation. *BMC Bioinformatics. 2016;17:73.

54. Barton NH. How does epistasis influence the response to selection? *Heredity (Edinb). 2016;118:96–109.

55. Stark C, Breitkreutz BJ, Reguly T, Boucher L, Breitkreutz A, Tyers M. *BioGRID: a general repository for interaction datasets. Nucleic Acids Res. 2006;34:D535–9.

56. Álvarez-Castro JM, Carlborg Ö. A unified model for functional and statistical epistasis and its application in quantitative trait loci analysis. *Genetics. 2007;176:1151–67.

57. Vitezica ZG, Legarra A, Toro MA, Varona L. Orthogonal estimates of variances for additive, dominance and epistatic effects in populations. *Genetics. 2017;206:1297–307.

58. Lee SH, van der Werf HJH. MTG2: an efficient algorithm for multivariate linear mixed model analysis based on genomic information. *Bioinformatics. 2016;32:1420–2.

59. VanRaden PM. Efficient methods to compute genomic predictions. *J Dairy Sci. 2008;91:4414–23.

60. Jaccard P. Nouvelles recherches sur la distribution florale. *Bull Soc Vaud Sci Nat. 1908;44:223–70.

61. Eynard SE, Windig JJ, Leroy G, van Binsbergen R, Calus MPL. The effect of rare alleles on estimated genomic relationships from whole genome sequence data. *BMC Genet. 2015;16:24.

62. Daetwyler HD, Capitan A, Pausch H, Stothard P, van Binsbergen R, Brøndum RF, et al. Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and complex traits in cattle. *Nat Genet. 2014;46:858–65.

63. Bolormaa S, Chamberlain AJ, Khansefid M, Stothard P, Swan AA, Mason B, et al. Accuracy of imputation to whole-genome sequence in sheep. *Genet Sel Evol. 2019;51:1.

64. Heidaritabar M, Calus MPL, Megens H-J, Vereijken A, Groenen MAM, Bastiaansen JWM. Accuracy of genomic prediction using imputed whole-genome sequence data in white layers. *J Anim Breed Genet. 2016;133:167–79.

65. Veroneze R, Lopes PS, Guimarães SEF, Silva FF, Lopes MS,HZarliusz B, et al. Linkage disequilibrium and haplotype block structure in six commercial pig lines. *J Anim Sci. 2013;91:3493–501.

66. Badke YM, Bates RO, Ernst CW, Schwab C, Steibel JP. Estimation of linkage disequilibrium in four US pig breeds. *BMC Genomics. 2012;13:24.
67. Andreescu C, Avendano S, Brown SR, Hassen A, Lamont SJ, Dekkers JCM. Linkage disequilibrium in related breeding lines of chickens. Genetics. 2007;177:2161-9.

68. Wientjes YCJ, Veerkamp RF, Calus MPL. The effect of linkage disequilibrium and family relationships on the reliability of genomic prediction. Genetics. 2013;193:621-31.

69. Uimari P, Tapio M. Extent of linkage disequilibrium and effective population size in Finnish Landrace and Finnish Yorkshire pig breeds. J Anim Sci. 2011;89:609-14.

70. Welsh CS, Stewart TS, Schwab C, Blackburn HD. Pedigree analysis of 5 swine breeds in the United States and the implications for genetic conservation. J Anim Sci. 2010;88:1610-8.

71. Vitezica ZG, Reverter A, Herring W, Legarra A. Dominance and epistatic genetic variances for litter size in pigs using genomic models. Genet Sel Evol. 2018;50:71.

72. Bulmer MG. The effect of selection on genetic variability. Am Nat. 1971;105:201-11.

73. Muir WM. Comparison of genomic and traditional BLUP-estimated breeding value accuracy and selection response under alternative trait and genomic parameters. J Anim Breed Genet. 2007;124:342-55.

74. Verrier E, Colleau JJ, Foulley JL. Long-term effects of selection based on the animal model BLUP in a finite population. Theor Appl Genet. 1993;87:446-54.

75. Wei M, Caballero A, Hill WG. Selection response in finite populations. Genetics. 1996;144:1961-74.

76. Mulder HA, Lee SH, Clark S, Hayes BJ, van der Werf JH. The impact of genomic and traditional selection on the contribution of mutational variance to long-term selection response and genetic variance. Genetics. 2019;213:361-78.

77. Liu H, Sørensen AC, Meuwissen THE, Berg P. Allele frequency changes due to hitch-hiking in genomic selection programs. Genet Sel Evol. 2014;46:8.

78. Fisher RA. The genetical theory of natural selection. Oxford: Oxford University Press; 1930.

79. Pszczola M, Strabel T, Mulder HA, Calus MPL. Reliability of direct genomic values for animals with different relationships within and to the reference population. J Dairy Sci. 2012;95:389-400.

80. Clark SA, Hickey JM, Daetwyler HD, van der Werf JHJ. The importance of information on relatives for the prediction of genomic breeding values and the implications for the makeup of reference data sets in livestock breeding schemes. Genet Sel Evol. 2012;44:4.

81. Lourenco DAL, Misztal I, Tsuruta S, Aguilar I, Lawlor TJ, Forni S, et al. Are evaluations on young genotyped animals benefiting from the past generations? J Dairy Sci. 2014;97:3930-42.

82. Weng Z, Wolc A, Shen X, Fernando RL, Dekkers JCM, Arango J, et al. Effects of number of training generations on genomic prediction for various traits in a layer chicken population. Genet Sel Evol. 2016;48:1-10.

83. Hill WG. "Conversion" of epistatic into additive genetic variance in finite populations and possible impact on long-term selection response. J Anim Breed Genet. 2017;134:196-201.

84. Huang W, Mackay TFC. The genetic architecture of quantitative traits cannot be inferred from variance component analysis. PLoS Genet. 2016;12: e1006421.

85. Cheverud JM, Routman EJ. Epistasis and its contribution to genetic variance components. Genetics. 1995;139:1455-61.

86. Manolio TA, Collins FS, Cox NJ, Goldstein DB, Hindorff LA, Hunter DJ, et al. Finding the missing heritability of complex diseases. Nature. 2009;461:747-53.

87. Marouli E, Graff M, Medina-Gomez C, Lo KS, Wood AR, Kjaer TR, et al. Rare and low-frequency coding variants alter human adult height. Nature. 2017;542:186-90.

88. Zeng J, de Vlaming R, Wu Y, Robinson MR, Lloyd-Jones LR, Yengo L, et al. Signatures of negative selection in the genetic architecture of human complex traits. Nat Genet. 2018;50:746-53.

89. Lehner B, Crombie C, Tischler J, Fortunato A, Fraser AG. Systematic mapping of genetic interactions in *Caenorhabditis elegans* identifies common modifiers of diverse signaling pathways. Nat Genet. 2006;38:896-903.

90. Tyler AL, Ji B, Gatti DM, Munger SC, Churchill GA, Svenson KL, et al. Epistatic networks jointly influence phenotypes related to metabolic disease and gene expression in diversity outbred mice. Genetics. 2017;206:621-39.

91. Strandén I, Christensen OF. Allele coding in genomic evaluation. Genet Sel Evol. 2011;43:25.

92. Wientjes YCJ, Bijma P, Vandenplas J, Calus MPL. Multi-population genomic relationships for estimating current genetic variances within and genetic correlations between populations. Genetics. 2017;207:503-15.

93. Wientjes YCJ, Calus MPL, Duenk P, Bijma P. Required properties for markers used to calculate unbiased estimates of the genetic correlation between populations. Genet Sel Evol. 2018;50:65.

94. Yang J, Benyamin B, McEvoy BP, Gordon S, Henders AK, Nyholt DR, et al. Common SNPs explain a large proportion of the heritability for human height. Nat Genet. 2010;42:565-9.

95. Vitezica ZG, Varona L, Legarra A. On the additive and dominant variance and covariance of individuals within the genomic selection scope. Genetics. 2013;195:1223-30.

96. Toro MA, Varona L. A note on mate allocation for dominance handling in genomic selection. Genet Sel Evol. 2010;42:33.

97. Su G, Christensen OF, Ostersen T, Henryon M, Lund MS. Estimating additive and non-additive genetic variances and predicting genetic merits using genome-wide dense single nucleotide polymorphism markers. PLoS One. 2012;7:e45293.

98. Sonesson AK, Wooll Williams JA, Meuwissen THE. Genomic selection requires genomic control of inbreeding. Genet Sel Evol. 2012;44:27.

99. Doekes HP, Veerkamp RF, Bijma P, Hiemstra SJ, Windig JJ. Trends in genome-wide and region-specific genetic diversity in the Dutch-Flemish Holstein-Friesian breeding program from 1986 to 2015. Genet Sel Evol. 2018;50:15.

100. Wooll Williams JA, Berg P, Dagnachew BS, Meuwissen THE. Genetic contributions and their optimization. J Anim Breed Genet. 2015;132:89-99.

101. Meuwissen THE. Maximizing the response of selection with a predefined rate of inbreeding. J Anim Sci. 1997;75:934-40.

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

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