Genomic prediction in the wild: A case study in Soay sheep

Genomic prediction, the technique whereby an individual’s genetic component of their phenotype is estimated from its genome, has revolutionised animal and plant breeding and medical genetics. However, despite being first introduced nearly two decades ago, it has hardly been adopted by the evolutionary genetics community studying wild organisms. Here, genomic prediction is performed on eight traits in a wild population of Soay sheep. The population has been the focus of a >30 year evolutionary ecology study and there is already considerable understanding of the genetic architecture of the focal Mendelian and quantitative traits. We show that the accuracy of genomic prediction is high for all traits, but especially those with loci of large effect segregating. Five different methods are compared, and the two methods that can accommodate zero-effect and large-effect loci in the same model tend to perform best. If the accuracy of genomic prediction is similar in other wild populations, then there is a real opportunity for pedigree-free molecular quantitative genetics research to be enabled in many more wild populations; currently the literature is dominated by studies that have required decades of field data collection to generate sufficiently deep pedigrees. Finally, some of the potential applications of genomic prediction in wild populations are discussed.


26
A major aim of evolutionary quantitative genetics is to measure and understand heritable genetic 27 variation, and to explain the maintenance of that variation in the face of natural and sexual selection 28 and genetic drift (Walsh and Lynch, 2018). One way that researchers have tried to answer these 29 questions is by conducting long-term ecological studies of natural populations (Charmantier et al., most (in) famous example of this is human height, where in a sample size of ~¼ million people, 697 59 genomewide significant loci only explained 16% of the heritability (Wood et al., 2014). Studies of natural 60 populations typically involve sample sizes of hundreds, or at most, a few thousand individuals, and 61 therefore statistical power to detect loci of medium-small effect is low. This has been recently 62 demonstrated in collared flycatchers, where association mapping with SNP chips and whole-genome 63 resequencing of individuals with extreme phenotypes failed to identify loci that explained variation in 64 the iconic sexually-selected trait, male forehead patch size (Kardos et al., 2016). The problem may be 65 exacerbated if phenotypes are hard to measure in the field or if environmental heterogeneity is a 66 further source of phenotypic variation. Put simply, even in very large experiments, the power to detect 67 most causal genes is very low. In natural populations, where sample sizes are currently smaller than in 68 studies of human or model organisms, the problem is far greater. This makes it impossible to describe 69 trait architectures or to understand the relationship between genotypes and phenotypes. How then, 70 can evolutionary quantitative geneticists use genomics tools to understand and predict the genetic 71 architecture and microevolution of their focal traits?

73
One possible solution to the low power / high polygenicity problem is to switch focus away from the 74 hunt for individual significant loci, and instead try to understand the genetic component of phenotypic 75 variation by using information from all of the typed SNPs. Two conceptually different approaches can 76 be taken. The first involves running an 'animal model' with the genetic relationship matrix estimated 77 from markers rather than a pedigree. It was recognised several decades ago that allele-sharing at 78 markers could be used to estimate relatedness between pairs of individuals, and that provided there 79 was variance in relatedness, the relatedness coefficients could then be regressed on phenotypic 80 similarity to predict heritabilities (Ritland, 1996, Ritland andRitland, 1996). Unfortunately, the reliable 81 estimation of relatedness, especially between distant relatives, requires several orders of magnitude 82 more markers than could realistically be typed at that point. As genotyping technology improved, 83 methods using marker-based relatedness became more tractable. In fact, marker-based approaches 84 are potentially more powerful than pedigree-based methods because the realised relatedness between 85 individuals rather than the expected relatedness can be estimated e.g. (Visscher et al., 2006, Hayes et 86 al., 2009b. Using a genomic relationship matrix (GRM) estimated from marker data, instead of the 87 pedigree-derived A matrix, it is possible to estimate quantitative genetic parameters in a traditional 88 mixed model Best Linear Unbiased Predictor (BLUP) framework (VanRaden, 2008). These models have 89 been used to estimate trait heritabilities in natural populations  90 2013, Bérénos et al., 2014, Silva et al., 2017 as well as in humans  and agriculturally 91 important organisms (VanRaden et al., 2009) . An extension of the model, where the GRM is estimated from a subset of markers, has been used to partition genetic variance into specific parts of the genome, distribution, allowing some markers to have large effects. In Bayes B, there are two distributions of 126 marker effects, with a proportion of markers assumed to have zero effect and the remainder with 127 effects drawn from a scaled-inverse χ 2 distribution. Following the development of Bayes A and Bayes B, 128 further refinements, such as the marker effects coming from different types of distribution, have been 129 made (Gianola, 2013, Habier et al., 2013, but the main principal of using a training population to 130 estimate marker effects, which are then used to predict phenotypes in a genotyped test population 131 applies.

133
The accuracy of genomic prediction can be estimated by correlation between GEBV and phenotype, 134 divided by the square root of the heritability (Meuwissen et al., 2013). Accuracy depends on a number 135 of factors such as the heritability (h 2 ) of the trait (Daetwyler et al., 2008, Goddard, 2009, Hayes et al., 136 2009a, the training population sample size (Habier et al., 2013, Meuwissen et al., 2001, the marker 137 density (Habier et al., 2009, Meuwissen and, the underlying genetic architecture (i.e.

140
Empirical investigation of these factors on genomic prediction accuracy have been described in the 141 animal and plant breeding literature (Heffner et al., 2011, Zhao et al., 2012, Hayes 142 et al., 2009a. Highly heritable traits generally have more accurate GEBVs because the environmental 143 and non-additive genetic effects are relatively small, meaning marker effects are more reliably 144 estimated. Increasing the size of the training population allows for more accurate estimate of marker 145 effects; consequently, prediction accuracy is also improved. Accuracy is generally greater with a higher 146 marker density, because markers will be in stronger linkage disequilibrium with unknown causal loci.

147
Conversely, accuracy tends to be lower in populations with larger effective population sizes, because 148 linkage disequilibrium between markers and causal loci is usually weaker. More formal descriptions of 149 how these factors affect accuracy are available elsewhere (Daetwyler et al., 2008, Goddard, 2009 150 Meuwissen et al., 2013, Hayes et al., 2009a.

152
The effects of trait genetic architecture and the method employed for the genomic prediction are 153 connected. Because the true genetic architecture, i.e. the actual distribution of effect sizes is generally 154 unknown, it is important that the posterior distribution of marker effects matches the true (but 155 unknown) distribution of SNP effects as closely as possible. Methods where all markers have an effect 156 drawn from a single normal distribution are unlikely to be good predictors of phenotypes determined 157 by a small number of large effect loci. Thus, methods that are flexible enough to capture different genetic architectures are potentially the most useful (Daetwyler et al., 2013). For example, methods 159 that model marker effects as a mixture of different distributions can accommodate large and small 160 effect loci (Erbe et al., 2012, Zhou et al., 2013.
While genomic prediction has become a widely used tool in animal and plant breeding, it has yet to be 163 widely adopted in quantitative evolutionary genetics studies of wild populations, even though there are 164 a number of long-term longitudinal studies where suitable phenotypic and genomic data have been 165 collected. Therefore, the main motivation of this study was to investigate the feasibility and accuracy 166 of genomic prediction in Soay sheep. Previous genetic studies of Soay sheep (Bérénos et al., 2015, 167 Gratten et al., 2007, Johnston et al., 2011 have demonstrated that some recorded 168 phenotypes have a simple Mendelian basis (e.g. coat colour, coat pattern), others are highly polygenic 169 (e.g. adult weight), and some are intermediate with many genes of small effect and a small number of 170 loci with quite large effect contributing to genetic variation (e.g. horn length). Thus, a second aim of 171 the study is to assess the accuracy of genomic prediction in different types of trait. Finally, the third aim 172 of the study is to compare the accuracy of different genomic prediction methods for traits with different architectures (Daetwyler et al., 2013, Meuwissen et al., 2013, de los Campos et al., 2013a. Thus, we

182
The Soay sheep is a primitive feral breed that resides on the St Kilda archipelago, off the NW of Scotland.

183
Since 1985, the population on the largest island, Hirta (57°48' N, 8°37' W), has been the subject of a 184 long-term individual based study (Clutton-Brock et al., 2004). Briefly, the majority of the sheep resident 185 in the village Bay area of Hirta are ear-tagged and weighed shortly after birth and followed throughout 186 their lifetime. Ear punches and blood samples suitable for DNA analysis are collected at tagging. Every

187
year during August, adult sheep and lambs are captured and morphological measurements are taken.

188
Winter mortality is monitored, with the peak of mortality occurring at the end of winter/early spring,  the genotype calling and quality control of the data are available elsewhere (Bérénos et al., 2015, 197 Bérénos et al., 2014, Johnston et al., 2016, Johnston et al., 2011. Briefly, pruning of SNPs and 198 individuals was performed using Plink v 1.9 (Chang et al., 2015). SNPs were retained if they had a minor

207
Phenotypic data 208 We studied phenotypic data collected for eight different traits: foreleg length, hindleg length, weight, 209 metacarpal length, jaw length, horn length (normal-horned males only), coat colour and coat pattern.
An aim of the study was to investigate whether comparisons between the alternative genomic prediction methods were sensitive to the underlying genetic architecture. Therefore, we chose to study 212 traits that had been the focus of previous gene mapping studies and were known to have variable 213 architectures (summarised in Table 1). Among the body size traits, foreleg length, hindleg length, weight and horn length were measured on live animals, whereas metacarpal length and jaw length 215 were taken from cleaned post-mortem skeletal material. Coat colour and coat pattern are independent 216 discrete traits that are recorded at the time of capture and remain fixed throughout life (Clutton-Brock 217 , Gratten et al., 2007. Further details of how the morphological traits 218 are measured can be found elsewhere (Bérénos et al., 2015, Bérénos et al., 2014, Johnston et al., 2011.

219
Analyses of live-capture morphological data were restricted to animals captured in August, that were 220 at least 28 months old, to remove most complications of growth and give consistency with previous 221 genetic studies of these traits (Bérénos et al., 2015, Bérénos et al., 2014. Male horn length included 222 measurements taken outside of August, as many males are captured and measured during the rut in 223 November. Post-mortem measurements were likewise restricted to animals who were at least 28 224 months old when they died.

226
Modelling non-genetic effects

227
The quantitative traits being investigated are known to vary between sexes and age classes and are 228 influenced by environmental effects (Bérénos et al., 2014). Because some of the software being tested 229 does not allow the inclusion of repeated measurements or fitting of random effects, it was necessary 230 to fit models with non-genetic effects prior to performing the genomic prediction. Therefore, mixed 231 effects models were run that adjusted phenotypes for fixed effects (sex and age of the animal) and nongenetic random effects (birth year, capture year and animal identity). Traits recorded on live animals 233 had repeated measurements, so individual identity was fitted to remove any permanent environmental 234 effect on phenotype. Traits recorded on skeletal material were recorded just once per animal. Models 235 were run using the lmer function in the R v3.6.0 package lme4 v1.1-23 (Bates et al., 2015). The random 236 effect of individual identity was extracted from each model and retained as the phenotype to be 237 analysed in the genomic prediction models. The two Mendelian traits are categorical traits that are 238 fixed throughout lifetime, unaffected by environmental conditions and with the same penetrance in each sex; therefore no phenotypic adjustments were necessary. Details of the model outputs are 240 summarised in Supplementary Table 1.

243
A major aim of this study was to compare different methods for obtaining genomic estimated breeding 244 values (GEBVs). Below we briefly discuss these methods, which rely on subtly different assumptions 245 about the distribution of marker effect sizes, and outline how each model was run.

247
The BayesA method was one of the two models in the first genomic prediction paper 248 (Meuwissen 2001). All markers have a non-zero effect size, and the marginal distribution of

264
This method of genomic prediction is similar to BayesA in that all SNP effects come from a single 265 distribution. However, BayesL models marker effects as following a double exponential 266 distribution rather than a t distribution. A consequence of using this distribution is that, relative   Table 1) and degrees of freedom set to 10. Note that the default parameters give very similar 293 values for genomic estimated breeding values, but we find that some samples of the MCMC 294 chain return zero or very low estimates of additive genetic and residual variance under the 295 default settings (see Supplementary Material Table S2). As with all analyses run in BGLR, the 296 BayesR analyses were run for a total of 120,000 iterations with a burnin of 20,000 and a 297 thinning interval of 10.

299
Because BayesR reports the genetic variance explained by each SNP, SNP effect sizes can be 300 aggregated (e.g. across the genome, across individual chromosomes) to describe the genetic 301 architecture of the trait. Detailed descriptions of the trait architectures are not the main goal 302 of this work, but we do report them in the Supplementary material. This is partially motivated by the knowledge that many of our focal traits already have well-described architectures (Table  are consistent with genome wide association studies. We also compare descriptions of trait 306 architectures between BayesR models using the default parameters and models using priors 307 informed by previous quantitative genetic studies of those traits (Table S2).

309
Assessing prediction accuracy and bias 310 Genomic prediction accuracy was assessed by cross-validation; each data set was randomly split into a 311 training population containing 95%, 50% or 10% of individuals and a test population containing the 312 remaining 5%, 50% or 90% of individuals. Twenty replicates were generated for every combination of 313 trait, model and training population size (8 traits * 5 methods * 3 training sizes * 20 replicates = 2400 314 models). For the six continuous traits, prediction accuracy was estimated as the correlation between 315 the GEBV and an individual's phenotype, divided by the square root of the heritability (i.e. rGEBV,y/h,

316
where h 2 is the heritability and y is the phenotype) Accuracy was averaged across the 20 replicates.

317
The phenotypes were the values pre-adjusted for age, sex, year of birth and year of death using the 318 same models as described in the section on Modelling non-genetic effects. Because the phenotypes 319 were adjusted for random and fixed effects, the heritability of each trait was estimated from the 320 adjusted phenotypes. These estimates of heritability will be greater than previously published 321 estimates, because non-genetic sources of variation have already been removed, but they are more 322 appropriate for measuring the observed and expected accuracy of the GEBVs. Trait heritabilities were 323 estimated using Animal models run in the R package MCMCglmm v2.29 (Hadfield, 2010). In these 324 models the additive genetic variance was estimated from a relationship matrix determined by the 325 population pedigree; the pedigree was constructed from a combination of behavioural observations 326 and genotype data from 315 SNPs. Pedigree construction methods are described elsewhere (Bérénos 327 et al., 2014). MCMCglmm models were run for 120,000 iterations with a burn-in of 20,000 and sampling 328 every 100 th iteration. In addition to recording the correlation between GEBVs and phenotypes, a linear 329 regression was performed to estimate the slope between the two variables. A slope of 1 is indicative of 330 there being no bias in GEBVs (Meuwissen et al., 2001, Moser et al., 2015. Slopes greater than one 331 would suggest that between-individual differences in GEBVs underestimate the difference in   (Table S2, Figure 1a). Unsurprisingly, the 350 accuracy declines when the training size is decreased. The accuracy of jaw length GEBVs were a little 351 lower than for the other traits. It is not immediately obvious why jaw length would have a lower accuracy, as it is not less heritable than the other traits, and although it is highly polygenic -no QTL 353 have previously been identified in GWAS studies -it is no different to adult weight in that regard 354 (Bérénos et al., 2015). The phenotypes of the two categorical traits were predicted with very high 355 accuracy, even when only 10% of animals were used in the training dataset (Table S2, Figure 1b). As 356 expected, for all traits, the accuracy was generally lower when the training set was smaller (Table S2, 357 Figure 1).

359
Comparisons between studies are not straightforward because accuracy depends on both biological 360 (e.g. genome size, effective population size, heritability) and technical (e.g. marker density, training 361 population size) factors (Daetwyler et al., 2010, Daetwyler et al., 2008, Goddard, 2009, de los Campos 362 et al., 2013b) , which inevitably will vary between studies. Certainly, the accuracy reported here is 363 comparable to situations in plant and animal breeding where genomic selection is routinely used for 364 trait improvement (Lin et al., 2014, Hayes et al., 2009a and is probably greater than has been observed 365 for morphological traits in commercial sheep populations (Auvray et al., 2014, Brito et al., 2017 . Therefore, the prospects for genomic prediction in this population are promising. If genomic prediction estimates of GEBVs are unbiased, then it is expected that regressions of GEBVs 370 on phenotypes should be equal to 1.0 (Meuwissen et al., 2001).There was a tendency for regression 371 coefficients to be a little less than 1.0, meaning GEBVs overestimated between-individual variation, but 372 this was largely driven by the models with smallest training sizes (Figure 2, Table S3). In the 50% and 373 95% training population datasets, there was relatively little bias with regression coefficients close to 1.0 374 for most traits and most models. Taking an average across 360 analyses (20 replicates x 3 training sizes 375 x 6 quantitative traits) BayesL had the mean regression coefficient closest to 1.0 (Table S3, top row).

376
However, of all of the methods, it had by far the largest standard deviation of regression coefficient, 377 and therefore may be the most prone to severe bias. Of the remaining methods, BayesR had the mean 378 regression coefficient closest to 1.0 and with the lowest standard deviation, although all of the other 379 methods performed quite similarly with respect to bias (Table S3).

382
Among most of the continuous traits, there was little difference in the accuracy (Table S2, Figure 1a).

383
This is not surprising as all of the models either make an assumption that the trait is polygenic, or can  (Moser et al., 2015, Erbe et al., 2012. From this point, we mostly focus on results 393 from the BayesR analyses, as that is the most accurate approach, and because it allows a more detailed 394 investigation of trait genetic architecture.

396
Comparisons between observed and expected accuracy

397
The relatively high accuracy of genomic prediction in Soay sheep is perhaps not especially surprising. It 398 is possible to predict the accuracy of genomic prediction if information about the effective population 399 size, genome size, recombination rate, trait heritability, and the genomic architecture are available 400 (Daetwyler et al., 2010. Here we used equation 2 of Daetwyler et al (2010): where rgĝ is the correlation between the true breeding value g and the GEBV ĝ. rgĝ is equivalent to the 403 rGEBV,y/h, the measure of accuracy we use in this paper (Meuwissen et al., 2013). Np is the number of  good fit between the observed and expected accuracy of genomic prediction (Fig 3); for the leg length 420 and weight traits, the observed accuracy was a little higher than expected.

421
The reasons for the accuracy of some traits being slightly better than expected are unclear. One possible 422 explanation is that the effective population size has been overestimated, which would result in a lower 423 predicted accuracy than if the correct Ne had been used. However, overestimating Ne would be 424 expected to affect the predicted accuracy of all traits, not just a subset of them. It is notable that 425 previous GWAS studies (Bérénos et al., 2015) of the three leg traits identified QTL of reasonably large effect (up to 15% of the additive genetic variance), despite the traits being reasonably polygenic. These 427 findings are supported by the BayesR analyses of genetic architecture (see Table S2). In general the 428 traits with larger effect size loci had the greatest accuracy, so perhaps the greater than expected 429 accuracy of leg length and weight traits is attributable, at least in part, to these relatively large effect documented (Gratten et al., 2007.

448
BayesR returns an estimate of the number of SNPs that explain the phenotypic variation. We urge 449 caution against taking these estimates as definitive (and note that the 95% confidence intervals span a 450 ~10-fold range for most traits; Table S4). However, the estimated number of SNPs was greatest for the 451 traits that were thought to be highly polygenic (weight, jaw length), and was least for the Mendelian

452
(coat colour and coat pattern) and major locus (horn length) quantitative traits. The estimated number

458
Chromosome partitioning plots ( Fig S1) were broadly consistent with previous work using chromosome-459 wide relatedness matrices (Fig 1 of Bérénos et al., 2015). However, in that study, two traits (adult weight 460 and jaw length) had significant correlations between chromosome length and the proportion of 461 variation explained by the chromosome, but in this study neither trait showed a significant relationship.

462
Nonetheless, the chromosomes that tended to explain the most variation in the earlier study, also 463 explained the most variation here. An exception was for jaw length, where Chromosome 20 explained 464 more than 30% of the additive genetic variance here, but did not in the earlier study. However, no SNP 465 on Chromosome 20 had a posterior inclusion probability > 0.5 ( Figure S2b), so there is no compelling evidence of a major locus QTL on that chromosome. Instead, posterior inclusion probabilities plots 467 indicate at least two regions of chromosome 20 contribute to jaw length variation (Fig S2b).

469
The in GEBVs between the two models was > 0.996 (Table S2). Genetic architectures were similar between 481 the two models, although there was a tendency for the heritability to be a little greater when using the 482 model that used informed priors (Table S2) In this population the accuracy of genomic prediction was comparable to that seen in applied animal 489 and plant breeding programs. In part, this is because Soay sheep have a history of isolation on a small 490 island and a relatively small effective population size. This means LD extends several megabases across 491 much of the genome (Bérénos et al., 2014, Feulner et al., 2013, and the number of SNPs (~38K) on a 492 medium density chip is sufficient to tag unknown causal loci. Future work will examine whether a higher

501
Instead, cross-sectional rather than vertical sampling of study populations is possible, which greatly 502 expands the range of organisms that could be studied. Genomic methods for describing genetic 503 architectures, like the BayesR approach, mean that many different aspects of genomic architecture (e.g. Hymenoscyphus fraxineus (Stocks et al., 2019). Being able to identify which young plants are most 515 resistant to dieback could be an essential tool in re-establishing populations devastated by disease.

516
Second, GEBVs can be used to better understand the 'invisible fraction' problem in evolutionary 517 biology, highlighted more than three decades ago (Grafen, 1988

525
Perhaps the most obvious application of genomic prediction in wild populations will be to explore 526 microevolutionary trends and in particular revisit the 'evolutionary stasis' problem (Merilä et al., 527 2001c). This refers to the observation that traits are frequently observed to be heritable and under 528 directional selection, yet they do not always evolve in response to that selection in the expected way.

529
Detailed descriptions and explanations of the possible explanations for stasis are available elsewhere (Kruuk et al., 2001, Merilä et al., 2001b, Merilä et al., 2001a, Merilä et al., 2001c. Quantitative genetic 531 approaches to understand evolutionary stasis have used EBVs derived from animal models applied to 532 pedigreed populations (Wilson et al., 2007, Coltman et al., 2003, Garant et al., 2004, but a number of problems and biases associated with this approach have been highlighted  the phenotype than true breeding values are (Postma, 2006). This is especially true for individuals 536 without many phenotyped relatives in the pedigree. This is problematic for studies looking at 537 microevolutionary trends in breeding values, as it means temporal changes in EBVs may reflect 538 environmental effects on phenotypes rather than genetic ones. Similarly, in a pedigree-based approach, 539 those individuals that lack phenotypic records and have few relatives have EBVs very close to the 540 population mean of true breeding values . In studies that explore temporal trends 541 in breeding values, these relatively uninformative individuals may be clustered towards either or both 542 ends of a time series. Both problems are largely overcome with GEBVs generated from a genomic 543 prediction test population, because each individual's genome rather than its phenotype or number of 544 phenotyped relatives is used to predict its breeding value. It has been shown through simulations that 545 genomic prediction is a more accurate method for estimating EBVs than pedigree-based methods when 546 the focal individual is not closely related to phenotyped individuals (Clark et al., 2012). Of course, other 547 problems highlighted in discussions of using breeding values to study microevolutionary trends, such 548 as failing to incorporate the uncertainty in GEBVs , must still be addressed, but 549 that is relatively straightforward if a Bayesian solution is used. Thus, it seems likely that genomic 550 prediction can pave the way for new analyses of microevolutionary trends.

552
In conclusion, this study shows that genomic prediction can reliably measure individual breeding values 553 in the Soay sheep population, and that high accuracy does not appear to be restricted to traits with 554 specific underlying genetic architectures. Approaches that can model zero effect, small effect and large 555 effect loci seem to be the most consistently reliable. Finally, we anticipate that similar studies will soon 556 be possible in many other previously under-studied organisms, paving the way towards both applied evolutionary quantitative genetics research and a re-exploration of some classic, yet unresolved, 558 problems. Agouti signalling protein (ASIP) --  N is the number of phenotyped and genotyped individuals for each trait. Genes that explain all or most of the variation in Oligogenic or Mendelian 562 traits and SNPs that explain more than 10% of the additive genetic variance of polygenic traits are all reported. References indicate publication 563 where estimation of heritability and linkage mapping and/or GWAS were performed. *Heritability estimates unique to this study are from a 564 MCMCglmm animal model using a pedigree-derived relationship matrix and phenotypes pre-adjusted for terms reported in Table S1. These 565 estimates are expected to be greater than the published estimates, due to the removal of some non-genetic sources of variation. However, they 566 are appropriate estimates when predicting the accuracy of GEBVs. 567 Figure 1: Accuracy of genomic prediction for (a) quantitative and (b) Mendelian traits. For quantitative traits, accuracy is measured as rGEBV,y/h, the mean correlation between GEBVs and the phenotype, divided by the square root of the heritability. For Mendelian traits, accuracy is determined as Area Under Curve from ROC analysis. All plots show the mean and SE obtained from 20 replicates of training sizes comprising 95%, 50% or 10% of the available individuals.