Two methods (Scheme A and Scheme B) were developed to optimize the relative weights on quantitative trait loci (QTL) and contributions of selected individuals simultaneously to maximize selection response while constraining the rate of inbreeding to the rate observed in gene assisted selection (GAS). In Scheme A, both the relative weights give to QTL and contributions of the selected individuals were optimized using a genetic algorithm. The possible solutions for relative weights of QTL and contributions of the selected individuals were encoded simultaneously. A physical selection population was used to evaluate the fitness of each encoded solution using stochastic simulation with 50 replicates. The fitness of each solution was the mean of all replicates for accumulative discounted sum of genetic means of all generations in physical selection population. In Scheme B, the optimization for relative weights on QTL was similar to Scheme A, and also was implemented based on a genetic algorithm. However, unlike Scheme A, an optimal contribution algorithm (OC) was used to optimize contributions of selection candidates. When compared with GAS, Schemes A and B resulted in up to 15.88 and 22.26% extra discounted sum of genetic value of all generations in a long planning horizon, respectively. Compared GAS+OC and Scheme B, most of the increase (about 78%) in genetic gain was produced by only optimizing contributions of selected individuals. The optimization for relative weight given to QTL just avoided the long-term loss (about 22%) observed in GAS scheme.

Introduction

With the rapid advancement of molecular genetics, many QTL/genes affecting quantitative traits in livestock species have been identified and mapped in recent years. Use of information on identified QTL/genes through gene assisted selection (GAS) in breeding program of livestock has been extensively studied (Dekkers, 2004, 2007). The common approach used in GAS was based on a simple index: I  =  α  +  EBV (Falconer and Mackay, 1996; Soller, 1978), where α is an estimated breeding value for the identified QTL and EBV is an estimated breeding value of the polygenic (excluding the marked QTL under evaluation) effects of an individual. In most of the published studies, GAS has implemented with equal emphasis on α and EBV (Abdel-Azim and Freeman, 2002; Gibson, 1994; Larzul et al., 1997; Pong-Wong and Woolliams, 1998; Ruane and Colleau, 1995) through truncation selection where all selected parents contribute equally to the next generation. It is generally concluded that extra responses are expected from GAS in the early generations of selection. However, these extra gains resulted from GAS cannot be maintained in a long run because the polygenic response that is lost in the initial generations could never entirely recovered in the later generations (Gibson, 1994). This phenomenon has been known as the Gibson effect.

Dekkers and van Arendonk (1998) developed a method to optimize selection on an identified QTL over a planning horizon of multiple generations, which solved the problem of lose of long-term response. Chakraborty et al. (2002) extended the method of Dekkers and van Arendonk (1998) to selection schemes where maximizing the sum weighted of short and long-term responses with different selection strategies for males and females, and multiple identified QTL allowing for non-additive (including dominance, epistasis and gametic imprinting) effects at the QTL. Tang and Li (2006a, b) extended the method of Chakraborty et al to allow optimization of selection on multiple traits with multiple QTLs in a population with overlapping generations. However, all of these previous studies assumed equal contributions of selected candidates with an infinite population sizes without consideration of the effect of accumulation of inbreeding.

Optimum contribution (OC) selection algorithms that maximize genetic response with a controlled inbreeding level have been developed during the last decade. Meuwissen (1997) and Grundy et al. (1998) suggested, independently, a method to maximize genetic gain while restricting the rate of inbreeding per generation to a predefined value in a population with discrete generations. The two methods were extended to populations with overlapping generations by Meuwissen and Sonesson (1998) and Grundy et al, respectively. Villanueva et al. (1999) used OC selection algorithms in an identified gene assisted BLUP scheme to maximize genetic gain with the constraint of inbreeding with equal emphasis to QTL and the polygene. However, these authors (Villanueva et al., 1999, 2002) also showed that the Gibson effect occurred in OC selection using a one-generation scheme (i.e., standard truncation selection and equal emphasis on the QTL and polygenic EBV in the selection criterion).

Mate selection is a breeding strategy that combines selection and mating simultaneously to maximize the benefit. Li et al. (2006, 2008) demonstrated the effectiveness of mate selection in exploiting non-additive QTL over multiple discrete generations in different breeding structures using a difference evolution method. This method was extended by Nishio et al. (2010) to a population with overlapping generation. Kinghorn (2010) proposed an algorithm to constrain mate selection efficiently, which extends the use of mate selection and enables implementation in relatively large programs across breeding units. However, this method does not consider polygenic effect in selection model, and neither controls the rate of inbreeding to a predefined level.

The strategy that optimizes the emphasis given to QTL and polygenic effects, and the strategy OC algorithm are implemented in different phases of breeding process. They can be combined together to further increase the selection response. Villanueva et al. (2004) tried to incorporate the two optimization methods in their gene assisted BLUP evaluation schemes, but their method was implemented through two steps. Firstly, the relative weight given to QTL and polygenic effects were optimized using deterministic simulation; secondly, OC selection was conducted based on previous results using the stochastic simulation. Under the two steps selection, the optimal result of the first step may not be suitable for the second step in a small population. Moreover, the two steps were not implemented simultaneously in a long planning horizon. The objective of this paper is to develop two methods that optimize the relative weight of QTL and contribution of individual simultaneously in a small population with a stochastic simulation and to further increases the selection response in a long-term planning horizon, while restricting the rate of inbreeding.

Materials and methods

Optimization objective

In general, a breeding objective is to maximize the sum of average genetic value of new progenies weighted by the corresponding weight in every generation within a planning horizon in each round of selection in multiple generations (Chakraborty et al., 2002):where G is the average genetic value for the t generation, w is the weight of the t generation , which was calculated based on the discount factor: w  = 1/(1 +  ρ), where ρ is the interest rate per generation (ρ  = 0.1, Li et al., 2006).

Optimization problem

In a common GAS strategy, there are two key problems that affect selection response of the terminal generation in a long planning horizon.

Gibson (1994) showed that common GAS strategy resulted in an extra selection response in earlier generations, but this extra response could not be maintained in the long-term selection due to the reason mentioned in the introduction section. The estimated genetic value for this strategy can be defined in the equation as follows:where q is the known genotypic value of the j QTL for individual i, is the estimated polygenic value (excluding the QTL under evaluation) of individual i. In this case, the weights are given equally to all of the QTL and to the polygene (i.e., all weights for QTL and polygene are one). This strategy can lead to the loss of polygenic response in the earlier generations of a long-term selection because the relative selection intensity given to QTL is larger than that to the polygene in early generations (Dekkers and van Arendonk, 1998). This implies that the weights on QTL and polygene are not optimal in the entire planning horizon.

Common GAS strategy utilizes the truncation selection to determine which individuals will become the parents of next generation, and the genetic contributions of all selected individuals of a specific sex were also assumed to be equal. However, in practice, individuals with different breeding value are expected to have different contributions at a certain rate of inbreeding. Thus, the genetic contributions of selected individuals in a conventional GAS strategy are not optimal. In addition, GAS strategy also pays no attention to the genetic relationship between individuals mated, and therefore, it couldn't control the average inbreeding level accumulated in the next generation.

Optimization strategy based on genetic algorithm

Genetic algorithms (Forrest, 1993) are a relatively new optimization technique which can be applied to the above optimization problems with multiple parameters. This technique does not ensure an optimal solution. However, it usually gives good approximations in a reasonable amount of time. Genetic algorithms are loosely based on natural evolution and use a “survival of the fittest” rule, where the best solutions survive and are varied until a desirable result is obtained. The basic process of genetic algorithm consists of encoding, evaluation, selection, crossover and mutation. In this study, for evaluating the fitness of every solution, a physical selection population over multiple generations was simulated randomly. On that basis, two different schemes (schemes A and B) were designed to optimize the previous two problems.

Simulation of physical selection population

In this study, the trait of interest was genetically controlled by an infinite number of additive loci (Bulmer, 1971), plus five identified quantitative trait loci (QTL). The genetic value of each individual was formulated as:where q denotes the j QTL genotypic value on individual i, u is the polygenic effect of individual i. The detailed parameters for simulating physical population are listed in Table 1 .

Table 1

The parameters of simulation.

Parameter	Value
Physical selection population
The phenotypic variance of interesting trait	1
The heritability of interesting trait	0.5
The proportion of genetic variation explained by all SNP markers	0.5
The number of simulated SNP marker	5
The initial favorable allele frequency of all SNP marker	0.2
The size of base and sub-generation population	500
The probability of sex	0.5
The selected proportion of sire	0.1
The selected proportion of dam	0.5
The number of simulated generation	20
The number of repetition	50
Encoding solution population
The size of encoding population	50
The number of encoded solutions to carry over unaltered to the next generation	5
The proportion of selected encoded solutions	0.5
The probability of mutation for a given encoded solution	0.3
The probability of crossover for selected encoded solution	1
The maximal iteration number	100

The base population was composed of 500 unrelated individuals (250 males and 250 females), and was simulated by one generation. The phenotypic variance (σ 2) was assumed to be equal to one. Here, for simplification, every QTL was supposed to explain equal additive genetic variance of the trait of interest. The initial favorable allele frequencies of all QTL were assumed as p  = 0.2. The genotype of each QTL of every individual was determined based on a uniform distribution U  ~ (0, 1)according to the initial favorable allele frequency p of QTL. The genotypic values of QTL were, and 0 for individual with favorable, unfavorable and heterozygous genotype, respectively (Falconer and Mackay, 1996), where σ 2 was the additive genetic variance of QTL. The polygenic effect u for each individual in base population was obtained from a normal distribution with a mean of zero and variance equal to h 2  − 5σ 2, where, h 2 is the heritability of trait under selection. The phenotypic value of an individual was calculated as: y  =  g  +  e , where e is a residual which was sampled from a normal distribution with a mean of zero and a variance of σ 2.

In subsequent generations, the offspring was generated by parents according to the mating system. The sex of each offspring was determined by U(0, 1) with equal probability of male and female. In this study, for simplification, all five QTL were assumed to be unlinked (i.e., they are independent each other). In this case, the first allele of all QTL for offspring was derived randomly from the two alleles of sire, and the second allele of all QTL was derived randomly from the two alleles of dam. The genotypic value of every QTL for each offspring was determined based on its genotype, and descended from the corresponding genotypic value in the base population. The polygenic effect of every offspring was generated based on the formula u  = 1/2u  + 1/2u  +  m , where u and u are the polygenic effects of sire and dam, respectively, and m is a Mendelian deviation which was sampled from a normal distribution with the mean zero and the variance of 0.5σ 2(1 − 0.5(F  +  F )), where σ 2 is the polygenic variance; F and F are the inbreeding coefficients of the sire and dam, respectively. The phenotypic value of every offspring was determined using the similar method in base population. The selection cycle was implemented over 20 generations.

Scheme A

In this strategy, the possible solutions for the relative weights given to all QTL and to the genetic contributions of all selected individuals were encoded simultaneously using the real number encoding (Fig. 1 ). Each encoded solution was a specific breeding strategy of previous simulated physical selection population. The encoded values of relative weight of QTL were applied at the evaluation phase of physical breeding process, and assigned a different weight to QTL and polygenic effects. The solutions of contributions of individuals were used to the mating system to allocate genetic contributions of selected individuals to next generation. A truncation selection method was used to determine the selected individuals of every generation in physical selection population, and the number of selected sires and dams were fixed in the whole breeding process.

Fig. 1

Scheme A based on genetic algorithm.

Encoding

Encoding relative weight of QTL (chromosome 1)

In previous simulated physical selection population, there were 5 QTLs, every QTL had 3 genotypes, and the round of selection was implemented 20 generations. Thus, a total of 300 (5 × 3 × 20) parameters for every solution were encoded. To maximize longer-term response, the emphasis given to QTL must be reduced because it reduces the selection pressure applied to polygene in the conventional GAS strategy (Dekkers and van Arendonk, 1998). Therefore, the boundary of every parameter was limited as [0, 1] in this study, which constrained the relative weight on QTL was smaller than the polygene in the optimization process. For every QTL, 3 random numbers were sampled randomly from the uniform distribution U(0, 1), and denoted the relative weight of favorable homozygote, heterozygous and unfavorable homozygote genotype, respectively.

Encoding genetic contribution of selection individuals (chromosome 2)

In previous simulated example, 25 sires and 125 dams were selected in every generation, and the cycle of selection was implemented 20 generations. Therefore, a total of 3000 (150 × 20) parameters were encoded, and the boundary of every parameter also was limited as [0, 1]. In addition, there is a constraint on the genetic contributions of physical selected individuals in a general mating system, i.e. the sum of genetic contributions of all selected individuals for specific sex was equal to one. Thus, for sires, 25 different random numbers were sampled firstly from a uniform distribution U(0, 1), and then these random number were divided by the sum of all random numbers () which caused the sum of all generated random number to equal one. For dams, same operation was implemented.

Evaluation

Like all optimization techniques, the genetic algorithm requires a measure of the utility of its proposed solutions. The fitness of every solution needs to be evaluated, and selection will be implemented according to the ranking order of fitness of all solutions (Forrest, 1993). In this study, a physical selection population of 20 generations (Section 2.3.1) was used to evaluate the fitness of every encoded solution (Fig. 1).

Assignment of relative weight of QTL

In Section 2.3.1, the genotype and polygenic effect of every individual had been simulated in the base generation. The selection index of every individual in the base generation was calculated as:where q is the known genotypic value of the j QTL based on specific genotype (favorable, heterozygous or unfavorable genotype),w is the relative weight of the j QTL which derived from the first 15 encoded values (the solution of first generation) of chromosome 1 for a specific encoded solution, is the estimated polygenic value, and calculated based on a BLUP mixed model:where  ' is the vector of corrected phenotypic value (y  ' =  y  −  q ), and are incidence matrices relating the observations to the fixed effects and polygenic effects, respectively, and is the vector of residuals. All individuals were selected based on their ranking order of selection index I according to the predefined selected proportion (10% for sire, 50% for dam).

Allocation of contributions of individuals and constraint of inbreeding

Genetic contributions of all selected individuals (25 sires and 125 dams) to next generation were determined based on the first 150 encoded values of chromosome 2 for specific encoded solution. The 150 encoded values were assigned to corresponding selected individuals based on a genetic algorithm, while restricting the rate of inbreeding. According to the description of Meuwissen (1997), for restricting the increase of inbreeding, the average coancestry between selected individuals should be controlled because the increase in coancestry equals the future increase in inbreeding. The average coancestry of selected individuals is equal to , where c is the vector of genetic contributions of the selected individuals to generation t  + 1, A is the matrix of additive genetic relationships among selected individuals in generation t, which equals twice the matrix of coefficients of coancestry. Thus, the objective function is to minimize the absolute value of difference between the average coancestry among selected individuals in generation t and the expected average inbreeding in generation t  + 1:where, is average inbreeding level among offspring in generation t  + 1, which equals to ΔF(t  + 1) , and ΔF is the desired rate of inbreeding which equals to the average rate of inbreeding observed in GAS. The absolute value is smaller, and denotes the effect that restricts the rate of inbreeding is better. The basic process of genetic algorithm was: 1) to generate an initial solution population with different permutation of c ; 2) to calculate the fitness of every solution based on the previous objective function; 3) to implement selection, crossover and mutation; 4) to generate new solution population. This cycle was repeated by 100 generations, and the best solution was used to assign the encoded contributions to the corresponding selected individuals in every generation.

Calculation of fitness

At last, the selected sires and dams were sampled based on their genetic contribution to generate a new sub-generation population. These previous breeding steps were repeated in the new sub-generation population based on the newly encoded value of the second generation for a specific encoded solution. The cycle of selection was implemented over 20 generations with 50 replicates for a specific solution (Fig. 1). The sum of the average genetic value of each generation weighted by corresponding weight in whole planning horizon was defined as the fitness of specific encoded solution (Section 2.1).

Selection

All encoded solutions were selected using standard truncation selection based on their fitness according to the predefined selection proportion (50%). In addition, there were 10% elite solutions that carried over to the next generation without crossover or mutation.

Crossover

Firstly, two mated solutions were determined whether to recombine based on a predefined probability of crossover (100%; Wright, 1991). Secondly, the positions of crossover were sampled from the uniform distribution U(1, 300) and U(1, 3000) for chromosome 1 and 2, respectively. Thirdly, a single point crossover was used to implement the recombination of encoded solutions. In this study, a SAS function GASETCRO (SAS/Genetics 9.2., 2008, IML procedure) was used to implement the crossover. The sum of new recombined chromosome 2 was re-constrained to one by using the value of every allele of chromosome2 divided by the sum of chromosome 2.

Mutation

In this study, a SAS function GASETMUT (SAS/Genetics 9.2., 2008, IML procedure) was used to implement the mutation. The function GASETMUT used a uniform mutation operator to change every value of chromosome 1 and 2 based on a varied rate of mutation and the predefined range of allele ([0, 1]), and then returned the compiled new chromosome 1 and 2. Here, the mutation rate was reduced progressively by 1% of initial rate of mutation (0.3; Wright, 1991) in every generation. The sum of new chromosomes 2 was re-constrained to one by the method described in Section 2.3.2.4.

Termination condition

Just as in nature, genetic algorithm operates the encoded population attempting to generate solutions with higher and higher fitness by repeatedly implementing the operations of crossover, mutation and selection until the termination condition is met. In this study, the termination condition was the difference value of the average fitness between previous and current generation was smaller than 0.0001 (Li et al., 2006).

Scheme B

In this scheme, only relative weights given to QTLs were encoded using the real number encoding (Fig. 2 ). The solutions of relative weight of QTL were used at the evaluation phase to assign corresponding weight to QTL and polygenic effect. The whole optimization process of genetic algorithm was similar to Scheme A except for the calculation of fitness. As Scheme A, a breeding process of 20 generations (Section 2.3.1) was used to evaluate the fitness of each solution. However, unlike Scheme A, an optimal contribution selection (Meuwissen, 1997) was used to obtain optimal number and genetic contributions of selected individuals based on a quadratic function in this strategy. This method can maximize the selection response of next generation while constraining average inbreeding coefficient to a predefined value. The optimal contributions of selected individuals in generation t were calculated as: where A - and EBV are the inverse of matrix of additive genetic relationship and vector of estimated breeding value in generation t, respectively, Q is a known incidence matrix for sex, λ and λ 0are Lagrange multipliers, respectively, is the average inbreeding among offspring in generation t  + 1, and calculated as , where ΔF is the desired rate of inbreeding which equals to the average rate of inbreeding observed in GAS. At last, the selected sires and dams were sampled based on their optimal contribution to generate the new offspring population. The loop of selection was implemented 20 generations with 50 replicates. The cumulative discounted sum of genetic value of all generations was defined as the fitness of specific solution.

Fig. 2

Scheme B based on genetic algorithm.

Comparison of selection strategy

When the termination condition of optimization scheme (Scheme A and Scheme B) based on the genetic algorithm was met, the encoded solution with the highest fitness in current generation is the optimal solution of optimization problem. For comparison, three types of selection schemes also were simulated, and they are described below: common gene assisted selection (GAS), GAS with optimal contribution (OC) control of inbreeding (GAS+OC) and conventional BLUP with OC control of inbreeding (CBLUP+OC).

GAS

In GAS, selection was directly on the QTL, and the effects of all QTL were known without error. The estimated breeding value (EBV) of every individual was calculated as: , where q is the known QTL genotypic value, is the estimated polygenic value, which was obtained from standard BLUP evaluation using the polygenic variance (σ 2) and the phenotypic values corrected for the QTL effect (Section 2.3.2.2). A standard truncation selection was used to select parent population, i.e., a fixed number of individuals (N males and N females) with the highest estimated breeding values were selected to be parents of next cycle. Each sire mated randomly to N /N dams.

GAS+OC

In this scheme, the estimation of EBV was the same with GAS. However, unlike GAS, an optimal contribution method (Meuwissen, 1997) was used to obtain optimal number and genetic contributions of selected individuals (Section 2.3.3) based on a constraint that the rate of inbreeding was equal to the average rate of inbreeding observed in GAS. After that, these selected sires and dams were sampled based on their optimal contribution to mate.

CBLUP+OC

Genetic evaluations were entirely based on the phenotypic and pedigree information. The EBV of individual was calculated from standard BLUP evaluation using the total additive genetic variance (σ 2  +  σ 2) and phenotypic values uncorrected. Like GAS+OC, a method described by Meuwissen (1997) was also used to optimize the contributions and number of selected individuals underlying the same rate of inbreeding observed in GAS (Section 2.3.3).

Results

The selection response and rate of inbreeding

The average selection responses of every generation for GAS, Schemes A and B are showed in Fig. 3 . As expected, both Schemes A and B resulted in more accumulative selection response than GAS at the same rate of inbreeding (Figs. 3 and 6 ). Compared the two components of selection response resulted from selection on QTL and polygene in Fig. 3, the accumulated selection responses from QTL selection improved more rapidly in early generations (the first nine generations) than that of polygenic selection in the GAS scheme. This phenomenon is similar to the results described by Gibson (1994). Contrarily, although Schemes A and B resulted in less QTL selection response than GAS in early generations, two strategies both maximized the early polygenic selection response by avoiding the loss of them. Moreover, they also maximized the QTL response (all QTL attain fixation) in the end of planning horizon (Figs. 3 and 4). The results of Schemes A and B were very similar with the results described by Chakraborty et al. (2002), which demonstrated the strategy that optimized the weight of QTL in this study substantially avoided the loss of early polygenic response. Compared Schemes A and B, the latter resulted more cumulative selection response than the former. In early generations, Schemes A and B obtained similar polygenic selection response, but the latter resulted more QTL response than the former. On the contrary, when QTL component of Schemes A and B attained fixation, the selection response from polygenic component in Scheme B improved more rapidly than Scheme A in late generations.

Fig. 3

Accumulative selection responses of QTL (﹍), polygene (﹎) and total (_) components for common gene assisted selection (○) and scheme A (●) and B (⋄).

Table 2 lists the genetic means of terminal generation and the cumulative discounted sum of all generations for Schemes A and B and GAS. The terminal genetic values of QTL component for three strategies were same because their favorable QTL alleles reached 1 in the end of planning horizon. The extra genetic values in the Schemes A and B were mainly resulted from polygene component, and the average polygenic genetic values in Schemes A and B were increased by 20.11 and 26.93%, respectively, at the terminal generation. Compared the results of the discounted sum of different response components in Table 2, the discounted sums of QTL in Schemes A and B were only accounted for 78.85 and 92.40% of the same discounted sum in GAS. These all due to the relative weight of QTL in the Schemes A and B were reduced in early generations and the polygenic discounted sum was improved by 33.78 and 36.69% in whole planning horizon and led to Schemes A and B resulted in 15.88 and 22.26% more discounted sum of overall response, respectively. Compared Schemes A and B, the latter resulted in more cumulative discounted response for QTL and polygenic component.

Table 2

The genetic value (mean ± SD) of terminal generation and cumulative discounted sum (mean ± SD) of genetic value of all generation for QTL, polygene and total components in common gene assisted selection (GAS) and Schemes A and B based on genetic algorithm.

Response component	Genetic value of terminal generation			Cumulative discounted sum of genetic value of all generations
	GAS	Scheme A	Scheme B	GAS	Scheme A	Scheme B
QTL	1.98 ± 0.00	1.98 ± 0.00 (99.79) a	1.98 ± 0.00 (100.00)	7.23 ± 0.41	5.70 ± 0.91 (78.85)	6.68 ± 0.45 (92.40)
Polygene	5.59 ± 0.21	6.71 ± 0.23 (120.11)	7.10 ± 0.25(126.93)	14.95 ± 0.92	20.00 ± 1.14 (133.78)	20.44 ± 0.93 (136.69)
Total	7.57 ± 0.21	8.69 ± 0.23 (114.79)	9.07 ± 0.25 (119.84)	22.18 ± 0.92	25.70 ± 1.37 (115.88)	27.12 ± 1.04 (122.26)

The elements in parenthesis are the ratio (%) of mean genetic values between Scheme A or B and GAS.

Table 3 also lists the total genetic means over multiple generations from GAS, GAS+OC, Schemes A and B, CBLUP+OC. In this study, for the purpose of comparison, GAS+OC, CBLUP+OC, Schemes A and B were implemented at the same rate of inbreeding observed in GAS (the rate of inbreeding in GAS was about 0.8%). As expected, optimization of genetic contribution of individuals obviously improved the genetic gain for all schemes. Compared GAS+OC and Scheme B, optimization for relative weight on QTL (OW) and contribution of selected individuals (OC) in Scheme B generated less response than only optimization for OC in GAS+OC in early generations (first seven generations). However, after 7th generation, Scheme B resulted in more and more selection response than GAS+OC because of implementing OW and OC simultaneously. The average extra responses for all generations resulted from only by implementing the OC (GAS+OC) and OW (the difference of genetic mean between Scheme B and GAS+OC) were 0.64 and 0.24, respectively. Most of the increase (78%) in gain resulted from the implementing of OC. The optimization for relative weight given to QTL just avoided the long-term loss of genetic variation (about 22%) observed in GAS. The terminal genetic mean of CBLUP+OC scheme was maximum, followed by Scheme B, GAS+OC and Scheme A. Contrarily, the sum of genetic mean of all generations of CBLUP+OC scheme was the least between them. The Scheme B resulted in maximum cumulative genetic mean for 20 generations, followed by GAS+OC and Scheme B.

Table 3

The total genetic means over multiple generations from gene assisted selection (GAS), GAS with optimal contribution (OC, Meuwissen, 1997) control of inbreeding (GAS+OC), Schemes A and B based on genetic algorithm, conventional BLUP with OC control of inbreeding (CBLUP+OC).

t	GASa	GAS+OCb	Scheme Bb	Scheme Ab	CBLUP+OCb
1	0.729	+ 0.249	+ 0.101	+ 0.028	− 0.129
2	1.417	+ 0.433	+ 0.184	+ 0.101	− 0.267
3	2.105	+ 0.672	+ 0.234	+ 0.135	− 0.375
4	2.773	+ 0.770	+ 0.431	+ 0.203	− 0.556
5	3.404	+ 0.684	+ 0.447	+ 0.257	− 0.706
6	3.986	+ 0.596	+ 0.556	+ 0.249	− 0.861
7	4.489	+ 0.398	+ 0.641	+ 0.366	− 0.819
8	4.932	+ 0.294	+ 0.704	+ 0.439	− 0.735
9	5.312	+ 0.314	+ 0.820	+ 0.520	− 0.583
10	5.664	+ 0.375	+ 0.902	+ 0.606	− 0.389
11	5.996	+ 0.445	+ 0.942	+ 0.685	− 0.234
12	6.312	+ 0.503	+ 1.040	+ 0.760	− 0.059
13	6.628	+ 0.521	+ 1.103	+ 0.797	+ 0.211
14	6.932	+ 0.652	+ 1.164	+ 0.866	+ 0.387
15	7.236	+ 0.771	+ 1.221	+ 0.918	+ 0.602
16	7.542	+ 0.926	+ 1.276	+ 0.933	+ 0.780
17	7.852	+ 0.978	+ 1.317	+ 0.950	+ 1.091
18	8.153	+ 1.004	+ 1.382	+ 0.955	+ 1.233
19	8.457	+ 1.003	+ 1.438	+ 0.962	+ 1.337
20	8.749	+ 1.106	+ 1.512	+ 0.965	+ 1.533
∑C	108.666	+ 12.695	+ 17.415	+ 11.694	+ 1.461

The GAS was firstly simulated. Then, GAS+OC, Scheme A, Scheme B and CBLUP+OC were constrained at the same rate of inbreeding observed in GAS, respectively.

The values of GAS+OC, Scheme A, Scheme B and CBLUP+OC were those deviated from GAS.

The last row was the sum of all generations.

The favorable allele frequency

Favorable allele frequencies of all QTL for GAS and Schemes A and B are showed in Figs. 4 and 5. Compared the GAS with the Scheme A in Fig. 4, the trend variation of favorable allele frequencies of the 5 QTL in GAS were very similar because they had same initial allele frequency in simulated schemes, and all favorable allele frequencies were increased rapidly to approach to fixation in first nine generations. However, the trends of the 5 QTL in Scheme A were fluctuant very differently and their allele frequencies were increased slower than that of GAS in early generations of the planning horizon. These results indicated Scheme A based on genetic algorithm reduced the selection pressure of the 5 QTL in early generations, and led all QTL to display slower fixation rates in the planning horizon. Compared Figs. 4 and 5, Scheme B caused favorable frequencies of some QTL (QTL2 and QTL3) to reach fixation more rapidly than GAS.

Fig. 4

Favorable allele frequencies of all QTL for common gene assisted selection (GAS, _) and Scheme A based on genetic algorithm (﹍).

Discussion

In this study, the main objective was to develop new optimization methods to maximize the sum of average genetic value of new progenies weighted by the corresponding weights in every generation of the planning horizon using a genetic algorithm. The results showed that two developed methods can result in substantial more selection response than common GAS when both the relative weight given to QTL and the genetic contribution of selected individuals were optimized simultaneously. The two methods also showed that the extra response could be maintained longer in the planning horizon.

Optimization of selection decisions in the original method of Villanueva et al. (2004) was implemented in two steps. Firstly, the relative weight given to QTL was optimized using the model described by Dekkers and van Arendonk (1998), which assumes equal contributions of selected individuals (i.e., truncation selection) and constant polygenic genetic variances in the whole planning horizon. Secondly, optimal genetic contributions of selected individuals were obtained based on the method of Villanueva et al. (1999). However, the first step of this method utilized deterministic simulation to describe the whole optimization process of relative weight on QTL, and the second step used stochastic simulation to optimize the contributions of selected individuals. Under such case, the optimal result of the first step may not be suitable for the second step in a small stochastic population. In this study, two different methods were developed to improve optimization strategies. In Scheme A, the relative weight of QTL and contributions of the selected individuals were optimized simultaneously based on a genetic algorithm. The relative weights of all QTL and contributions of selected individuals were first encoded simultaneously (chromosome 1 and 2). Then, they were changed simultaneously by crossover and mutation operators in the optimization process of genetic algorithm. Therefore, they were implemented based on one step in the random optimizing process. In Scheme B, like Scheme A, the optimization of relative weight of QTL was also implemented based on a genetic algorithm. However, different from the Scheme A, an optimal contribution algorithm described by Meuwissen (1997) was used to control the rate of inbreeding in Scheme B, and then determine the selected numbers of sires and dams and the optimal contributions of selected sires and dams (Figs. 1 and 2). In this method, the contributions of selected individuals did not need to be encoded in advance, and the numbers of selected sires and dams were varied. The major difference of Scheme A is that the numbers of selected sires and dams in Scheme A were fixed because of the fixed selected proportion of sires and dams and fixed population size, more importantly, also because the contributions of selected individuals needed to be encoded together with relative weight of QTL in advance for each solution in genetic algorithm. Therefore, the Scheme A maximized the cumulative discounted sum of genetic value of all generations subjected to the constraint that the numbers of selected sires and dams were fixed in every generation, and it was a local optimization strategy relative to Scheme B. In other words, the ability that utilized the optimization of contribution of selected individuals to increase the gain in Scheme A was not better than Scheme B. This reason may be used to explain why Scheme B resulted in more selection response than Scheme A (Fig. 3, Tables 2 and 3).

In this study, two different strategies (Schemes A and B) were developed to optimize the relative weight given to QTL and contributions of selected individuals simultaneously at the same rate of inbreeding in GAS. For constraining the inbreeding, Scheme A utilized the process of allocation of contributions of selected individuals to restrict the increase of inbreeding. The contributions of selected individuals were encoded in advance. A genetic algorithm was used to minimize the absolute value of difference between the average coancestry among selected individuals in generation t and the expected average inbreeding in generation t+1 by permuting the order of encoded solution of genetic contribution of selected individuals, i.e., found an optimal order to assign the encoded contributions to corresponding selected individuals to restrict the average coancestry among selected individuals to an expected level (which control the increase of inbreeding of next generation) in every generation. The object for restricting the inbreeding in Scheme A was the average coancestry of selected sires and dams. In Scheme B, the method for constraining the inbreeding was also based on average coancestry information using a quadratic optimization (Meuwissen, 1997). However, unlike Scheme A, the object for restricting the inbreeding in Scheme B was the average coancestry among all selection candidates, moreover, the contributions of all candidates were not determined in advance. Thus, Scheme B had a larger search space of contributions of individuals than Scheme A. This maybe is another reason why Scheme B resulted in more selection response than Scheme A (Fig. 3, Tables 2 and 3).

Villanueva et al. (2004) defined “long-term” as the generation where the favorable alleles of all QTL were fixed in all schemes compared. In this study, the relative weight of the 5 QTL were limited as the range of [0, 1] in the optimal selection based on the genetic algorithm. Here, a base precondition is that the favorable alleles of these 5 QTL will reach fixation before the end of a planning horizon in GAS (the favorable alleles of all 5 QTL for GAS reached fixation in first 9 generations in this study) because relative weight given to QTL ([0, 1]) in optimal selection will delay their paces towards to the fixation of the favorable QTL alleles (Fig. 2). The objective of this optimal selection based on genetic algorithm was to control the increasing paces of frequency of favorable allele of all 5 QTL in the early generations and to reach their fixation at the last generation of the planned horizon (Figs. 3 and 4). This optimization process is tried to avoid the loss of polygene selection responses component in the earlier generation meanwhile maximized the total selection response in the planning horizon. However, long-term selection will likely to increase the number of QTL under selection gradually, under such situation, the favorable QTL alleles will likely not be fixed in limited generations (such as 20 generations). The optimal selection may improve the relative weight of QTL to maximize the discounted sum of genetic value of all generations in a planning horizon. The range of relative weight of QTL also may exceed the [0, 1] in the optimization process. In general, the occurrence of the long-term loss depends mainly on the proportion of genetic variation explained by polygene in GAS. Smaller proportion, less long-term loss will be happened in the terminal generation of planning horizon (Dekkers et al., 2002; Dekkers and Chakraborty, 2001). When all genetic variations of interesting trait are explained by markers (such as genomic selection), theoretically, there are no loss of response in a very long-term. However, the generation number that all favorable alleles attain fixation will be very long under this case. In practical breeding program, the general breeding scheme is to maximize the accumulative discounted response in a planning horizon with limited number of generation. Therefore, when going for short accumulative gain, little emphasis will be given to these QTL with low frequency and small effect, which will cause some loss of response in current genomic selection scheme (Sonesson et al., 2010). The method developed in this paper based on genetic algorithm can optimize many parameters simultaneously, and may be used to optimize relative weights of many markers in genomic selection. Of course, in such situation, the calculation cost will be very very expensive.