Many authors (Goddard, 2009; Jannink, 2010; Sonesson et al., 2012; Clark et al., 2013; Sun et al., 2013) have expressed the importance of reducing inbreeding in PS and GS for long-term success. They argued that GS is likely to lead to a more rapid decline in the selection response unless new alleles are continuously added to the calculation of GEBVs, stressing the importance of balancing short and long term gains by controlling inbreeding in selection.

Let A be a matrix of pedigree based numerator relationships or the additive genetic relationships between the individuals in the genetic pool (this matrix can be obtained from a pedigree of genome-wide markers for the individuals) and let c be the vector of proportional contributions of individuals to the next generation under a random mating scheme. The average inbreeding and co-ancestry for a given choice of c can be defined as r=12c′Ac. If b is the vector of GEBV's, i.e., the vector of BLUP estimated BV's of the candidates for selection. The expected gain is defined as g = c′b. Without loss of generality, we will assume that the breeder's long term goal is to increase the value of g.

In Wray and Goddard (1994), Brisbane and Gibson (1995), and Meuwissen (1997) an approach that seeks minimizing the average inbreeding and co-ancestry while restricting the genetic gain is proposed. The optimization problem can be stated as

where ρ is the desired level of gain.

This problem is easily recognized as a Quadratic optimization problem (QP). There are many efficient algorithms that solves QP's so there is in practice little difficulty in calculating the optimal solution for any particular data set. Recently, several allocation strategies were tested using QP's in (Goddard, 2009; Schierenbeck et al., 2011; Pryce et al., 2012). It is easy to extend these formulations to introduce additional constraints as positiveness, minimum-maximum for proportions, minimum-maximum for number of lines (cardinality constraints).

Some authors recommended mate selection approaches that also seek a balance between gain and inbreeding from an animal breeding perspective (Allaire, 1980; Jansen and Wilton, 1985; Kinghorn, 1998). Simple methods such as sequential selection (Pryce et al., 2012) or linear programming (LP) (Jansen and Wilton, 1985; Weigel and Lin, 2000) have been used to find mate designs to avoid more-related pairs and find less-related pairs. Some approaches example the use of genome-wide marker information (Fernández et al., 2001; Pryce et al., 2012; Sun et al., 2013), but these approaches mainly deal with the restricted case of mating after evaluation of breeding population and selection. Kinghorn, in a series of articles (Kinghorn, 1998; Kinghorn and Shepherd, 1999; Kinghorn, 2011), describes an algorithmic approach that separates the optimization and the objective function for the mate selection approach and therefore can be used for a wide array of optimization criteria (MSI) with hard and soft constraints. The MSI recommended by Kinghorn (1998) involves calculation of the merit of potential calves as half of the sire's breeding value plus half the dam's breeding value minus a weighted penalty on the estimated progeny inbreeding coefficient. The method in Pryce et al. (2012) uses a similar MSI based on phased haplotypes. These sources also introduce a number of constraints related to constraints on reproduction or resources. The optimization criteria in these papers have the same flavor as the quadratic optimization problem in Equation (1) and is expected to give similar solutions with the parental contributions approach for large breeding populations. We agree with these previously recommended approaches in the sense that mate selection is an optimization problem. However, we use markers, genomic prediction models and estimates of these models from phenotypic experiments (as the training data becomes available) and include a genomically estimated additive variance term for each mate pair; therefore the optimization problem in the recent paper is a related but different one; for instance, our criterion is neither linear nor quadratic.

By solving the QP in Equation (1) for varying values of ρ, or using the similar criteria in the mate selection approaches, we can trace out an efficient frontier curve, a smooth non-decreasing curve that gives the best possible trade-off of genetic variance against gain. This curve represents the set of optimal allocations and it is called the efficiency frontier (EF) curve in finance (Markowitz, 1952) and breeding literature.

The implementation of PS in our simulations did not use any genotypic information or pedigrees. Basically, it referred to selecting the individuals with best observed phenotypes to be parents in the next generation. Results elsewhere (Forni et al., 2011) indicate that there would be no significant differences between PS and GS if a pedigree from many generations is used in pedigree based estimation of BV's. In addition, there are methods to combine pedigrees with marker based relationship matrices (Legarra et al., 2009; Meuwissen et al., 2011) which would result in a yet another selection approach.

Based on the selection intensity the best individuals were determined in each step based on their phenotypic values and a random mating was employed among these individuals to determine the mates. On the other hand, GS used GEBVs obtained from models that are updated by phenotypic values every two cycle using only the most current genotypic and phenotypic data. The selection of parents in each cycle was done according to the GEBVs and the selection intensity. The mate assignment was at random as in PS. Efficient GS referred to the selection of parents in each cycle using the optimal parental contribution proportions obtained by solving the optimization problem in Equation (1). The specific solution in each cycle was selected among the Pareto optimal solutions along the EF curve by setting ρ to the 90th percentile of the GEBVs in the current population (we have decreased ρ if a gradually if the solution was dominated by a single parent). The mates for next cycle were assigned using the optimal contribution proportions.

There are several alternative measures of inbreeding based on mating plans (Leutenegger et al., 2003; Wang, 2011). In this article, we have used a measure derived under the infinitesimal genetic effects assumption proposed by Quaas (1988) and Legarra et al. (2009). A measure of gain, i.e., the total expected breeding value of the progeny, can also be calculated from the results of the same authors (Quaas, 1988; Legarra et al., 2009). However, in our belief, the expected value by itself is not a good measure of possible gains since it carries no information about the variability of BV's among full-sibs. Therefore, we have derived a measure called the risk of a mating plan (this is related to the concept of “usefulness”) by increasing the expected of the progenies by a small amount (the intensity is controlled by the parameter λ₁) proportional to their expected variance (standard deviation) calculated under the infinitesimal effects assumption.

Combining the measures of inbreeding and risk into one leads to the formulation of the mating problem:

where λ₂ ≥ 0 is the parameter whose magnitude controls the amount of co-ancestry in the progeny, and the minimization is over the space of the mating matrices P₃₂ construction of which is described in detail below. λ₁ controls allele heterozygosity weighted by the marker effects and λ₂ controls allele diversity. When λ₁ = 0 the risk measure is the same as total expected gain.

Now, we give the details of how the measures Risk(λ₁) and Inbreeding are defined in this paper. First, we make the following assumptions: (a) Diploid behavior at meiosis, (b) Uncorrelated genes distribution, (c) Absence of allelic interactions, (d) No multiple alleles at those loci controlling the character, (e) No of genotype-environment interaction.

Let b=(b1′,b2′,b3′)′ denote the vector of genetic effects corresponding to the parents and progeny, where b₁ and b₂ are the genetic effects of the N parents and b₃ are the genetic effects of the N_c progeny. Let the pedigree based numerator relationship matrix for the individuals in b be A and this matrix is partitioned as

corresponding to the partitions of b. Suppose, we also have the markers for the parents in the second partition, and u₂ = Ma where M is the matrix of minor allele dosages, coded as 0, 1, and 2. Let M_c be the N × m marker allele frequency centered incidence matrix and a be the vector of marker effects. Variance-covariance of b₂ can be written as

with G, the genomic relationship matrix, defined as McMc′k where k=∑j=1m2pj(1-pj) is twice the sum of heterozygosities of the markers (VanRaden, 2008).

Following Quaas (1988) and Legarra et al. (2009), let P be a matrix containing the transitions from ancestors to offspring. We will refer P as the mating or parentage matrix. Then, we can write b = Pb + ψ where ψ is the vector of Mendellian samplings and founder effects with a diagonal variance D. In particular, using only the rows of P corresponding to the b₃ the relationship is written as

which can also be stated as a regression equation of the form b3=(I-P33)-1(P31b1+P32b2+ψ3) (Quaas, 1988). The variance-covariance matrix of b₃ is given by

The variances caused by Mendelian sampling in D₃ are related to inbreeding in the parents via

where F₁ and F₂ are the inbreeding coefficients of the two parents which can be extracted from the diagonals of G. The variance-covariance formula reduces to

if all the founders are genotyped (no P₃₁), and a relatively simple mating strategy is assumed where founders are the only parents and no back-crossing is allowed (P₃₃ = 0). This is the assumption made for the remainder of this paper and in this case P₃₂ is a N_c × N matrix (N_c children from N parents) with each row having two 1/2 values at positions corresponding to two distinct parents or only a value of 1 at the position corresponding to the selfed parent. All the other elements of this matrix are zero. Nevertheless, one can easily imagine situations where some of the founders are not genotyped or when some of the progeny also have progeny, then the formula in Equation (3) will be relevant. If some founders are not genotyped but a pedigree is available relating them to the rest of the founders then the variance-covariance for the founders, Var(b₁, b₂), can be calculated using the relationship matrix in Legarra et al. (2009). Furthermore, construction of the mating matrices for more complex mating plans is described in Quaas (1988).

Var(b₃) gives us the expected variance-covariance of the progeny given the mating matrix P₃₂ and the realized relationship matrix G of the parents. This can be used as to measure the expected genetic diversity of a mating plan: we can use a measure in line with the inbreeding term in Equation (1) by

We also need a measure for genetic gain. A simple measure of gain for a given mating plan expressed in P₃₂ can be constructed from the expected value of b₃:

and an overall measure can be written as

Finally, we want to complement the measure “gain” with a measure of within cross-variance for the genetic levels of children of the parent pairs. Suppose the organism under study is diploid. We can re-code the markers matrix M coded as −1, 0, and 1 into a N × m matrix M^* using the information in the marker effects vector a such that markers are coded as the number of beneficial alleles as 0, 1, or 2. This is achieved by first obtaining marker effects estimates and then using the sign of the estimates to determine what is a beneficial allele. We can also obtain a related marker effects vector a^* by replacing the original marker effects by the effects of the beneficial alleles (a^* = |a|) so that we have Ma=(M*-1N×m)a*. For a given parent pair, we can calculate the vector expected number of beneficial alleles of the children of these parents using a transition vector p as μ = E(m) = p′M^*. In addition, for each locus we can calculate the variance for the number of beneficial alleles from the number of alleles the parents have and put them in a vector which we will denote by σ_p = (σ_p1, σ_p2, …, σ_pm). Calculation of elements of σ_p from the coding in M^* can be as in Table 1. We define risk measure for this parent pair as

where λ₁ ≥ 0 is the risk parameter and m is the number of markers. The risk of a mating plan (which is expressed in P₃₂) is the sum of all the risk scores for all mate pairs in that plan which we will denote by Risk(P₃₂, λ₁).

If the risk parameter λ₁ is set to zero then we have

The magnitude of λ₁ is related to the desire of the breeder to take advantage of within cross variances and encourages mating parents that are heterozygotes at quantitative trait loci (QTL).

In this sense, the efficient mating problem can be stated as an optimization problem as follows:

In the above optimization problem, we are trying to minimize the inbreeding in the progeny while the risk is set at the level ρ ≥ 0. In the remainder of this paper, we will use the the equivalent formulation of the mating problem in Equation (2).

The risk measure for each mate pair is calculated as the sum of the estimated breeding value for the progeny plus λ₁ times the expected standard deviation of the BV's of the progenies of this pair. When λ₁ is set to zero, only the expected BV's for the progenies are taken into account. If λ₁ > 0, then the variation expected in the progenies enters the picture and the pairs that have higher expected variance among their progenies will have higher risk values among mate pairs that have same expected breeding value. Keeping everything else constant the effect of an increased variance for the progenies of a mate pair is increasing the chances that one of the progenies of this mate pair will have a high breeding value compared to the mean value of the progenies. For instance, if we assume a normal distribution for the breeding value of the progenies of a mate pair, about 16% of the progeny will have breeding value higher than one standard deviation above the mean, i.e., mate pairs that have high variance have the potential to produce progenies with high BV's. However, since the distribution is symmetric around the mean, the same mate pair will have equal chances of producing progenies with low BV's. Therefore, we named our measure as “risk”. We note also that other measures of expected variance have been proposed (Zhong and Jannink, 2007) and these can be used instead of the one we have proposed above. For example, it is possible to calculate this variance by simulating progenies and their expected BV's for parent pairs, and one can easily include information about the LD in these simulations. Our choice of the measures of inbreeding and risk are driven by computational efficiency in order to keep the optimization over the mates feasible. In the next section, we give some examples of GM; we provide a procedural form of the GM algorithm used in these examples by Procedure 1.

The differences between the expressions in Equation (1) (which can be written in the current context by replacing A with G as c′Gc) and the current criterion in Equation (2) should be quite noticeable: c is the vector of parental contribution proportions and P₃₂ is the mating matrix that leads to progeny and this reflects the shift of focus to selection of mates rather than providing only proportions but leaving the mating problem unanswered. Given P₃₂ proportion of parental contributions can be calculated, however there is no way of obtaining an optimal P₂₃ from the knowledge of c. This focus shift also leads to the introduction of the diagonal matrix D₃ for the variances in progeny caused by Mendelian sampling of the alleles in the parents in the latter. In addition, a representation of the criterion in Kinghorn (Kinghorn, 1998) using the current notations is given by MSI=1Nc′P23GP23′1Nc+λ1Nc′P23b. Writing c^′=1Nc′P23 and MSI=c^′Gc^+λc^′b, we can see that this and Equation (1) have the same form and since we can expect these criterion to give similar results if N_c is large and no other constraints were imposed.

As compared to the quadratic optimization problem in Equation 1 which admits easy solutions, the optimization problem in Equation (2) is a combinatorial problem whose order increases with the number of individuals in the breeding population and the number of progeny. The list of all possible mates including self-mates is a list of length N_p(N_p + 1)/2. The problem can be stated as selecting N_c elements with replication from this list. The number of un-ordered N_c-tuples of an N_p(N_p + 1)/2-set is (Np(Np+1)/2)Nc/Nc!. We have devised a genetic algorithm to find good solutions for this optimization problem in reasonable computing time. The optimization procedure that is used in optimal mating problem is a modified genetic algorithm (GA) supplemented with tabu search (See Procedure 2).

Genetic algorithms (Holland, 1973; Davis, 1991; Goldberg, 2006) are particularly suitable for optimization of combinatorial problems. The idea is to use a population of candidate solutions that is evolved toward better solutions. At each iteration of the algorithm, a fitness function is used to evaluate and select the elite individuals and subsequently the next population is formed from the elites by genetically motivated operations like crossover and mutation. Tabu search is a search where most recently visited solutions are avoided by keeping a track of the previously tried solutions. This avoids many function evaluations and decreases the number of iterations till convergence, it is especially useful for generating new solutions toward convergence.

Each point in the solution space of the optimal mating problem can described by a frequency distribution over the the set of all possible mate pairs. In our case, the recombination of two solutions in this GA are obtained by summing the frequency distributions of these two solutions and sampling with probabilities corresponding to this combined frequency distribution. With a given probability of mutation, a mutation event follows the crossover event and is achieved by decreasing the frequency of a mate that has positive frequency by some integer value less than the current frequency of that mate and increasing the frequency of some other mate pair by the same amount. In our algorithm, the average magnitude of change that is allowed during mutation is called mutation intensity.

It should be noted that the solutions obtained by any GA or any other evolutionary algorithm may be sub-optimal and different solutions can be obtained given a different starting population. Another layer of safety is obtained if the algorithm is started from multiple initial populations and an island model of evolution is used where separate populations are evolved independently for several steps and then the best solutions from these algorithms becomes the initial solutions to evolutionary algorithm, this strategy is also partially employed.

We included the minimum number of parents as a parameter: “minparents” in long term simulations using GM. This allowed us to run the simulations many times without interference. Nevertheless, a better approach in practical situations would be to plot the whole frontier surface and select a solution that has a good risk to diversity ratio.

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