The Modular Breeding Program Simulator (MoBPS) allows efficient simulation of complex breeding programs

Keywords: resource management, animal breeding, Göttingen Minipigs, optimum contribution selection, OGC, OCS.

Breeding programs aim at improving the genetic characteristics of livestock populations with respect to productivity, fitness and adaptation, while controlling negative effects such as inbreeding or health and welfare issues. As breeding is affected by a variety of interdependent factors, the analysis of the effect of certain breeding actions and the optimisation of a breeding program are highly complex tasks. Although quantitative theory, such as the breeders’ equation (Falconer and Mackay 1996), can give approximations when considering single generations or relatively simple designs, this typically is not sufficient for todays’ complex breeding programs.

With increasing computing power, a relatively new approach to answer these questions is the execution of stochastic simulation by tools such as MoBPS (Pook et al. 2020). Stochastic simulations typically simulate each individual mating, including recombination events in the meiosis, and common breeding actions such as the phenotyping on a per individual basis. In comparison to a deterministic analysis based on quantitative genetic theory, such as implemented in ZPLAN+ (Täubert et al. 2010), stochastic simulations are computationally highly demanding, especially with respect to the simulation of meiosis and the estimation of breeding values.

While simulation tools have become more and more sophisticated and powerful, both research questions of scientists and practical design questions of breeders have become more complex at the same pace. Thus, simulation software frameworks constantly need to grow to allow for a larger population, to provide more modelling options for genetic traits, and to include new selection techniques, while still being accessible to a wide audience and not overwhelm the potential user by the sheer flood of parameter options. To help with those issues, Simianer et al. (2021) recently proposed a unifying concept to describe breeding programs in a comprehensive, unambiguous and reproducible way.

Most large-scale breeding programs use genotyping of at least some individuals to get a more accurate breeding value estimation, to control inbreeding rates, to confirm the pedigree, and to potentially reduce the generation interval (Meuwissen et al. 2001; Schaeffer 2006; Henryon et al. 2014). However, in small-scale breeding programs, genotyping still represents a substantial cost factor – and often it is not immediately clear whether the cost of genotyping is compensated by the benefits of genomic analysis.

In the present work, we will introduce various new features in the simulation framework MoBPS to more realistically model modern breeding programs and will display the potential of the framework to answer real-world breeding questions through simulation studies. To this end, we will exemplarily analyse the potential benefits of genotyping in the breeding program of Göttingen Minipigs (Reimer et al. 2020). Göttingen Minipig is the smallest commercially available minipig breed under a controlled breeding scheme, its breeding is globally organised in multiple isolated and independent colonies, and it is used in biomedical research as an animal model for various human disease areas and in preclinical toxicology and safety assessment studies (Pedersen and Mikkelsen 2019). The interested reader is referred to Simianer and Köhn (2010) for more details. Furthermore, we will introduce a new graphics processing unit (GPU)-based implementation to derive the genomic relationship matrix and to perform genomic best linear unbiased prediction (GBLUP) by solving the mixed model with known variance components, which massively reduces the computing time of this procedure.

In this study, the R-package MoBPS (Pook et al. 2020) was used as the backend simulation software. To generate a realistic initial population, real genomic data of 150 Göttingen Minipigs genotyped with an Affymetrix Axiom™ Porcine Genotyping Array (600K) was used as a set of founders. On the basis of these individuals, a stock of an arbitrarily chosen size of 5000 animals was generated by randomly mating founder individuals. Furthermore, three generations of random mating within the flock were simulated to build up a basic pedigree structure. Founder individuals were assumed to be related according to a genomic relationship matrix derived following VanRaden (2008). In case no genomic data are available, a valid alternative would be the simulation of a far higher number of generations, say 1000, to build up a basic LD structure (function: founder.simulation() in MoBPS). For all matings, an average litter size of 5.6 was assumed, with the distribution of the litter size given in Table 1. For simplicity reasons, litter size was assumed to be independent from genetic components.

Table 1.  Distribution of the simulated litter size

Three independent quantitative traits were simulated assuming 1000 underlying purely additive quantitative trait loci each, with effect sizes drawn from a Gaussian distribution and heritabilities of h₁² = 0.3, h₂² = 0.3 and h₃² = 0.1 respectively. Traits were standardised to have an underlying true genetic variance of 100 and all animals were phenotyped for all three traits. Trait 1 can be thought of as a main growth trait (such as ‘smallness’, since a small body size is highly desirable for a porcine animal model), and Trait 3 represents a typical functional trait such as fertility. Phenotypes of the second trait were assumed to be observed in a discrete distribution to mimic a trait such as the teat count, with the distribution in the founder population given in Table 2. Note that the underlying trait still has a quantitative background and only the phenotypic expression is affected by this. An in-depth description of the generation of discrete value traits can be found in the MoBPS user guidelines (available at https://github.com/tpook92/MoBPS). In case real genomic and phenotypic data would be available, one could consider estimating quantitative trait loci effects via ridge regression BLUP to place marker effects onto known-effect regions (effect.estimate.add() in MoBPS).

Table 2.  Distribution of the teat count in the founder population

A schematic overview of the breeding program itself is given in Fig. 1. The given breeding scheme can also be found as a template in the MoBPSweb application at www.mobps.de (Pook et al. 2021). In the following, we will consider five different selection strategies that require different amounts of genotyping:

Fig. 1.  Schematic overview of the baseline breeding program. This figure is a visualisation generated via the MoBPSweb graphical interface (www.mobps.de; Pook et al. 2021).

1. Phenotypic selection

2. Pedigree BLUP (no genotyping)

3. Single-step genomic BLUP (10% genotyping)

4. Single-step genomic BLUP (25% genotyping)

5. Genomic BLUP (100% genotyping)

In case a breeding value estimation was performed, the animals of the last two generations were used with a pedigree-depth of 7. The genomic relationship matrix for the genotyped individuals was derived according to VanRaden (2008), and in cases in which not all individuals were genotyped, the combined relationship matrix for genotyped and non-phenotyped animals was derived via the single-step framework (Aguilar et al. 2010; Legarra et al. 2014).

The top 50 breeding boars and top 500 breeding sows were selected on the basis of a selection index, with equal weights on the three quantitative traits.

As traits were modelled to be uncorrelated and have the same overall genetic variance, single-trait models were used for the breeding value estimation and the selection index was a plain combination of the estimated breeding values. To handle more complex and realistic scenarios, both modelling of correlated traits and more complex selection indices that utilise the reliability of the estimated breeding values and correlation between traits are implemented in MoBPS (Hazel and Lush 1942; Miesenberger 1997).

Furthermore, we considered the use of optimum contribution selection (Meuwissen 1997) for the selection of the breeding boars. For this, the method introduced in Wellmann (2019), which is implemented in the R-package optiSel (Wellmann 2017), was used with the target function to maximise the breeding value (‘max.BV’), while limiting the increase in average kinship (‘ub.sKin’). This functionality is also implemented within the MoBPS framework and was tested as an alternative selection criterion in all scenarios, except plain phenotypic selection.

The described breeding cycle was repeated for 10 generations and averages of 50 independent simulations per scenario were used as results.

To reduce computing times, variance components were assumed to be known in the mixed model. Furthermore, the miraculix package was extended to utilise GPU resources. This GPU version relies on the cutlass and cuSolver libraries by NVIDIA (see https://github.com/NVIDIA/cutlass) for optimising memory latency and utilises bit compression to minimise memory transfers. By combining tailored two-bit matrix multiplication algorithms with instruction sets of state-of-the-art GPUs further optimises the actual computing times required for the breeding value estimation.

Across the different scenarios, we observed increasing prediction accuracy when a more sophisticated model was used for genomic prediction and/or the share of genotyped animals was increased. For example, the relative gain after 10 breeding cycles was 60.0% higher for the main growth trait when the phenotypic selection was replaced by a pedigree-based breeding value estimation (Fig. 2a). The relative gain when switching from pedigree BLUP (Scenario 2) to GBLUP (Scenario 5) was still substantial (15.3%), while further gains when using single-step BLUP were low. The genetic gain for the teat count trait was slightly lower than for the purely quantitative trait with the same heritability (Fig. 2b). Note that the discrete distribution was simulated with a maximum value of 18, which is reached by a substantial share of the animals (~10%, Fig. 3), hindering differentiation between top animals and, thus, also further genetic progress. Details on the obtained prediction accuracies are given in Fig. 4. As the additive genetic variance is going down slightly over time, while residual variances are simulated to remain the same, prediction accuracies slightly decrease over time. In practice, the introduction of new genetic variance and a larger training population over time will usually compensate for this.

Fig. 2.  Average increase in underlying true genomic value for the sow cohort over time for (a) the main growth trait, (b) teat count, and (c) the fertility trait. 95% confidence bands are given for the average value and based on 50 independent runs each.

Fig. 3.  Distribution of the teat count traits for Scenario 2 at different time points for the sows cohort: (a) Year 1, (b) Year 5 and (c) Year 10.

Fig. 4.  Accuracy of the breeding value estimation of time for (a) main growth trait, (b) teat count, and (c) fertility trait. 95% confidence bands are given for the average value and on the basis of 50 independent runs each.

In regard to analysing genetic diversity, average levels of kinship and inbreeding (Fig. 5) were calculated on the basis of the share of the genome in identity-by-descent (Donnelly 1983). As underlying points of recombination are known in MoBPS, this can directly be calculated. Note that this measure does not account for inbreeding occurring before the founder generation of the simulation and thus should be seen as an indicator for the increase of inbreeding only since the founder generation. Since the overall selection pressure is lowest when applying phenotypic selection, this results in by far the lowest increase in both kinship and inbreeding. Although genetic gains were highest in the scenario using genomic BLUP (Scenario 5), increases in inbreeding are lower than for all other BLUP-based scenarios, with GBLUP leading to a reduction of the average kinship of 24.8%, 22.5% and 18.3% compared with Scenarios 2, 3 and 4 respectively. Since GBLUP, by design, should be superior to pedigree-based methods to distinguish between close relatives (e.g. full sibs) and usually leads to the selection of animals from more families (Daetwyler et al. 2007), this should not be that surprising.

Fig. 5.  Development in average (a) kinship and (b) inbreeding for the sow cohort over time. 95% confidence bands are given for the average value and on the basis of 50 independent runs each.

The performance of the optimum contribution selection for selecting breeding boars highly depends on the model parameter ‘ub.sKin’, which determines the maximum increase in average kinship per generation and can be interpreted as a weighting between genetic progress and genetic diversity. Across all scenarios, results were very similar; therefore, only Scenario 5 (GBLUP) is presented here. Overall, there are weightings that lead to both higher genetic gain and lower average kinship, hence being strictly better than the plain selection of the top breeding boars. For example, with similar genetic progress, we have a reduction in the average kinship by 15.2%, from 0.101 (plain selection) to 0.086 after 10 generations (Fig. 6). When restricting the increase in kinship to 0.005 per generation, this results in a genetic gain of 43.3 and an average kinship of 0.091 (compared with 42.3 and 0.101 in Scenario 5). When allowing for an extremely high increase in kinship, the short-term genetic gain is slightly increased, but after 10 generations, resulting gains are lower as genetic diversity has decreased heavily.

Fig. 6.  The realised average kinship versus the genetic gain after 10 generations of optimum contribution selection (line) against Scenario 5 (cross). Realisation that are both right and below the non-OGS scenario (cross) indicate strict improvements of the breeding programs in regard to kinship and genetic gain (red area).

Depending on the scenario, a single simulation run took up to 21.9 min on a single core of an Intel Xeon Gold 6154 3.00GHz CPU. Thereof, 10.6 min were used for the generation of new animals, calculation of pedigree, single-step and genomic relationship matrices took 1.3, 3.2 and 0.7 min, and solving of the mixed model equations took 3.8 min. The whole process of a single breeding value estimation for a single trait including the calculation of the relationship matrix for 10 000 full genotyped individuals took 30 s (or 14 and 9 s when using 5 and 20 cores). The GPU implementation on a RTX 2080 Ti reduced the 30 s to less than 4 s (Schlather et al. 2021). In comparison, rrBLUP (Endelman 2011), which estimates also the variance components, would need 198 s. For a larger number of individuals, the traditional rrBLUP approach took ~100 times as long as the GPU-based approach via miraculix (Table 3). Although almost negligible in our scenarios, the computing times for solving the mixed model scales cubically with the number of individuals. Therefore, this is a critical bottleneck for large-scale simulation, as other computational factors such as the generation of individuals, scale approximately linear. The assumption of known heritability might not be true in practice; however, it is common practice to avoid re-estimating variance components for every new breeding value estimation. Especially for large-scale breeding programs, heritability estimates should be very accurate, so that the impact on the prediction accuracies should be only minimal.

Table 3.  Computing time in seconds for a single breeding value estimation, including the calculation of the genomic relationship matrix and solving of the mixed model (as performed in Scenario 5), depending on the number of included animals and the implementation

As expected, the use of genomic data improves the accuracy of genomic prediction and thus leads to higher genetic gain in the breeding program of Göttingen Minipigs. In a simulation study, we could quantify the expected gain to be ~15% between a pedigree-based and a fully genomic evaluation, while reducing gains in inbreeding by 25%. Note that one of the main breeding objectives in the Göttingen Minipig is a uniform reaction in medical tests and, thus, high genetic similarity among animals is actually desired, while negative effects of inbreeding must still be avoided. This is a particularly difficult breeding objective, as low genetic diversity usually means high inbreeding and vice versa.

The benefits of genotyping go even further, as it allows to validate the pedigree data and enables further downstream analyses such as genome-wide association studies. Furthermore, genotyping should be of even higher benefit when estimating breeding values for non-phenotyped animals to potentially save phenotyping costs and/or reduce the generation interval (Schaeffer 2006). Assuming genotyping costs of 25€ per animal, the yearly costs of genotyping will sum up to 125 000€. A breeding company needs to offset these costs against the potential benefits. Note that MoBPS could also answer follow-up questions and/or be used to assess more complex genotyping strategies by using different marker density arrays, use of (pooled) sequencing, and similar, which were neglected here. Furthermore, we can confirm that the use of optimal contribution selection can lead to strict improvements over the plain selection of top breeding boars with equal contributions and thus should be implemented in any case.

As a more general note, the present study has provided an exemplary case on the use of stochastic simulation to evaluate and compare breeding strategies. Without any doubt, our simulation study is a simplification of reality, as the underlying traits are usually much more complex, both in terms of genetic architecture and in how they are derived (e.g. the weight trait in the Göttingen Minipig is a weighted combination of multiple weight measurements at different age points; Köhn et al. 2007). As this, in the end, should still result in a trait that approximately corresponds to a Gaussian distribution, this should only have a small impact on the overall results.

In general, before setting up a simulation study, the most important factors influencing the simulation should be identified. The inclusion of any such factor should be well considered, as it might unnecessarily increase the complexity of the model. Note that many model extensions will affect different scenarios in similar ways. Thus, the estimation of the overall effects might be slightly off; however, the relative differences among scenarios will typically still be accurately estimated. As this is usually of interest in a simulation study, we generally recommend keeping a simulation study on an appropriate level of simplification. However, if a more complex design is needed after all, modern simulation software such as MoBPS provides a variety of features for modelling, such breeding programs (Büttgen et al. 2020) in a flexible but still highly efficient way.

The new features introduced to the MoBPS simulation framework allow for both more realistic modelling of common breeding actions (e.g. optimum contribution selection, adaptable litter sizes) and trait architectures (e.g. discretely distributed phenotypes), which, to our knowledge, are not implemented in any other stochastic simulation framework. Note that a variety of other features have been added since the initial release of MoBPS (Pook et al. 2020) to handle subpopulations, calculate average offspring performance, use different arrays for genotyping or features to automatically avoid fullsib/halfsib matings that are not relevant for the present study but further enhance the design options in the MoBPS framework. In combination with the graphical user interface (Pook et al. 2021), the variety of pre-implemented methods allows the simple use of the MoBPS framework not only for research but also for users without scripting experience, in industry, or to assist teaching of breeding program design. At the same time, the flexibility and modular structure allows researchers to assess new and experimental methodology, such as new selection techniques or a new breeding value estimation approach. MoBPS then only takes care of the background processes (e.g. meiosis, trait simulation, data storage) and enables the simulation of an even wider range of breeding programs. Additionally, the high computation efficiency of MoBPS in regard to both highly efficient storage (Pook et al. 2020) and faster computing than with other simulators enables the efficient simulation of a high number of runs of various scenarios.

Just a single simulation run for each considered scenario will usually not be enough to detect differences between designs unless they are of large effect. In contrast to real-world breeding programs, it is, however, possible to repeat simulations without requiring a large amount of time or financial effort (besides extended computational needs). Furthermore, external factors that potentially affect the outcome do not need to be controlled. This enables the extended statistical analysis of minor changes of a breeding scheme. In contrast to most deterministic simulations (Täubert et al. 2010) or analyses based on quantitative genetics and selection theory (Falconer and Mackay 1996), the use of stochastic simulation further allows for the estimation of the variability of the outcome. Variability can be accounted for by the calculation of confidence intervals, but also the direct analysis of individual simulations to make sure a certain breeding target is reached with a given probability. In particular, for large-scale breeding programs, mimicking all computational steps performed in practice for each generation in each run of each scenario computationally can be very costly, thus highlighting the need for both simplifications and highly efficient software. For the breeding value estimation, we recommend the use of known heritabilities, as resulting prediction accuracies are only mildly affected and variance components are at least approximately known (and in many practical breeding programs are not even estimated in each cycle either). As computing time of the breeding value estimation scales cubically, further simplifications, such as the use of approximate algorithms, such asthe algorithm for proven and young (Misztal 2016), or the down-scaling of cohort sizes, might be considered.

The script used to perform the simulations can be found at https://github.com/tpook92/MoBPS/tree/master/Exemplary_scripts. Genomic data used in the simulations is available on reasonable request to the corresponding author (Torsten.pook@uni-goettingen.de).

The authors declare no conflicts of interest.

This research did not receive any specific funding.