About Us | Help Videos | Contact Us | Subscriptions

The Plant Genome - Original Research

Accuracy of Genomic Prediction in a Commercial Perennial Ryegrass Breeding Program


This article in TPG

  1. Vol. 9 No. 3
    unlockOPEN ACCESS
    Received: Nov 01, 2015
    Accepted: May 13, 2016
    Published: August 4, 2016

    * Corresponding author(s): dfe@dlf.com
Request Permissions

  1. Dario Fè *ab,
  2. Bilal H. Ashrafb,
  3. Morten G. Pedersena,
  4. Luc Janssb,
  5. Stephen Byrnec,
  6. Niels Roulunda,
  7. Ingo Lenka,
  8. Thomas Didiona,
  9. Torben Aspc,
  10. Christian S. Jensena and
  11. Just Jensenb
  1. a DLF A/S, Research Division, Højerupvej 31, 4660, Store Heddinge, Denmark
    b Aarhus Univ., Dep. of Molecular Biology and Genetics, Center for Quantitative Genetics and Genomics, Blichers Allé 20, 8830 Tjele Denmark
    c Aarhus Univ., Dep. of Molecular Biology and Genetics, Crop Genetics and Biotechnology, Forsøgsvej 1, 4200, Slagelse, Denmark
Core Ideas:
  • High accuracies for genomic prediction in a perennial ryegrass breeding program
  • The additive genetic variance can be traced by genotyping assays
  • Predictions work across different generations and in different traits
  • Good prospects for the implementation of genomic selection in perennial ryegrass


The implementation of genomic selection (GS) in plant breeding, so far, has been mainly evaluated in crops farmed as homogeneous varieties, and the results have been generally positive. Fewer results are available for species, such as forage grasses, that are grown as heterogenous families (developed from multiparent crosses) in which the control of the genetic variation is far more complex. Here we test the potential for implementing GS in the breeding of perennial ryegrass (Lolium perenne L.) using empirical data from a commercial forage breeding program. Biparental F2 and multiparental synthetic (SYN2) families of diploid perennial ryegrass were genotyped using genotyping-by-sequencing, and phenotypes for five different traits were analyzed. Genotypes were expressed as family allele frequencies, and phenotypes were recorded as family means. Different models for genomic prediction were compared by using practically relevant cross-validation strategies. All traits showed a highly significant level of genetic variance, which could be traced using the genotyping assay. While there was significant genotype × environment (G × E) interaction for some traits, accuracies were high among F2 families and between biparental F2 and multiparental SYN2 families. We have demonstrated that the implementation of GS in grass breeding is now possible and presents an opportunity to make significant gains for various traits.


    CRR, crown rust resistance; FR, fructan; G × E, genotype × environment; GBLUP, genomic best linear unbiased prediction; GBS, genotyping-by-sequencing; GEBV, genomic estimated breeding values; GS, genomic selection; LD, linkage disequilibrium; MAS, marker-assisted selection; NDF, neutral detergent fiber; PP, parental population; QTL, quantitative trait loci; SNP, single-nucleotide polymorphism; SY, seed yield; SYN2, multiparental synthetic; TKW, 1000-kernel weight; YLST, trial within year, location, and score; YLT, trial within location and year

Perennial ryegrass is the most cultivated forage species in temperate regions across northwest Europe, America, South Africa, Japan, Australia, and New Zealand (Humphreys et al., 2010). It is an obligate allogamous species (Cornish et al., 1979) bred in genetically heterogeneous families. Breeding programs for this crop started around the 1920s, and since then, breeders have mainly relied on phenotypic recurrent selection of superior individuals (Wilkins and Humphreys, 2003; Conaghan and Casler, 2011; Hayes et al., 2013). Through decades, this approach has led to significant improvements in several characters such as rust resistance, spring growth, and aftermath heading. However, poor gains have been achieved in traits like dry matter production or seed yield (Sampoux et al., 2011). Furthermore, current breeding programs tend to be expensive and time consuming, needing from 3 to 5 yr for a selection cycle and up to 15 yr for the release of a new cultivar (Casler and Brummer, 2008; Posselt, 2010).

New perspectives were opened in the last decade by the implementation of marker-assisted selection (MAS). Analyses performed on perennial ryegrass revealed the presence of several quantitative trait loci (QTL), uncovering associations between different traits and certain DNA regions (reviewed in Shinozuka et al., 2012). Marker-assisted selection approaches were used in many crop species, having a positive impact on breeding for qualitative traits. However, this technique was less effective in improving complex traits (Moreau et al., 2004; Bernardo, 2008; Xu and Crouch, 2008), as such traits are generally affected by many genes each with a small effect. In many cases, those effects are too small to be detected. That results both in a large number of false negatives and in the occurrence of the so called Beavis effect (Beavis, 1998; Xu, 2003), which consists in a bias of the allelic effect estimation because the estimated QTL effects are actually sampled from a truncated distribution.

These problems may be overcome by the use of GS (Meuwissen et al., 2001), which select the best individuals and families based on their genomic estimated breeding values (GEBVs). The calculation of GEBVs is performed by using all markers at the same time, and it is known as genomic prediction. Implementation of GS into breeding programs is now possible through development of new technologies that allow the deployment of high-throughput genotyping assays in a cost-effective manner. Using a high marker density, a large part of the causal polymorphisms is expected to be in linkage disequilibrium (LD) with at least one marker. In GS, the LD can be also tracked across families, enabling the computation of marker effects estimates at population level, while traditional MAS only allows such estimates at family level (Jannink et al., 2010). Furthermore, using all markers means that GS has the potential to explain the whole, or a very large part of, total genetic variance, thus overcoming the Beavis effect and leading to unbiased estimates of the breeding value (Hayes et al., 2013).

The practical application of genomic prediction implies different steps. In the first step, statistical models are developed and tested on a training set of individuals or families that are both genotyped and phenotyped. Such a population should contain as much variation as possible, minimizing the level of relationship between individuals (Pszczola et al., 2012). In the following step, GEBVs for other individuals or families (validation set) are estimated based on their genomic relationships with the training set and without the need to phenotype them.

Genomic selection has become a commonly used technique in animal breeding (e.g., Hayes et al., 2009). In crops, it has been widely investigated on simulated data (reviewed in Jannink et al., 2010). The first studies on empirical data were published in 2009 for maize (Zea mays L.) and barley (Hordeum vulgare L.) (Lorenzana and Bernardo, 2009; Piepho, 2009). Later studies also included wheat (Triticum aestivum L.) and other species (reviewed in Lin et al., 2014). The main focuses of these studies were accuracy and comparison between different statistical methodologies. Overall, GEBV predictions for plants were significantly superior to the ones based solely on phenotype. That was true both for simulated and empirical data (Jannink et al., 2010; Lin et al., 2014) especially in the case of traits with low heritabilities and with a large training dataset. Furthermore, models worked well for both polygenic and oligogenic traits (Iwata and Jannink, 2011), also in strategies based on multitrait indexes (Heffner et al., 2009). Reduction in phenotyping is expected to significantly decrease the length of the breeding cycles, resulting in higher genetic gains per year (Jannink et al., 2010) and, potentially, in a significant drop in costs of developing new lines. Heffner et al. (2010) estimated the expected annual gain to exceed present gains by about three times in maize and two times in winter wheat. However, there are still concerns about the impact of an extensive use of GS over the long term (Heffner et al., 2009; Jannink et al., 2010; Jannink, 2010; Nakaya and Isobe, 2012), which, at the present state of technology, make it difficult to completely replace the traditional breeding strategies with GS.

Regarding forage species, the potential to implement GS in ryegrass has been recently explored by Hayes et al. (2013). Other studies have also investigated the potential for GS in model species, white clover (Trifolium repens L.), tall fescue (Festuca arundinacea Schreb.), and Phalaris spp. (Simeão Resende et al., 2014; Forster et al., 2014). Hayes et al. (2013) define perennial ryegrass as an optimal target species for introducing GS because of the polygenic nature of many of its commercial traits, which are usually measured in the later part of the breeding cycle. However, they also emphasize the presence of three relevant problems: (i) nonextensive LD, (ii) likely high effective population size because of the outbreeding nature of the species, and (iii) need of radical changes in the current breeding schemes. The overall positive expectations were confirmed by Fè et al. (2015a) who showed the superiority of genomic prediction over MAS using heading date as model trait with a very high heritability. However, the literature still lacks extensive studies performed on a wider set of traits characterized by lower heritabilities.

The aim of this study is to investigate the possibility for implementing GS in perennial ryegrass breeding by using data from a commercial breeding program. The breeding material consisted of families produced by DLF A/S. Traits included in this report are: (i) resistance to crown rust (caused by Puccinia coronata Corda f.sp. lolii Brown), which is, so far, the most serious foliar disease affecting ryegrasses; (ii) seed yield; (iii) 1000-kernel weight; and (iv) quality parameters measured as fiber and fructan content.

Materials and Methods

Plant Material, Genomic, and Phenotypic Data

Data were collected from a total of 1918 families of perennial ryegrass, derived from the breeding program at DLF A/S (Store Heddinge, Denmark). The plant material consisted of two different sets of families:

Set 1 included 1791 biparental F2 families produced between 2000 and 2012 with parents chosen from 198 parental populations (PPs). The PPs were selected from a seed bank that included breeding material from DLF as well as varieties available on the market produced by DLF or other companies. The breeding procedure involved pair crosses between single plants from different PPs (each single plant was used only in one pair cross) followed by hand-harvesting from both parent plants, pooling and multiplication of the F1 seeds, and harvesting of F2 seeds as described in Fè et al. (2015b).

Set 2 included 127 multiparental synthetic (SYN2) families, each obtained by crossing five to 11 single plants from different F2 families originating from different PPs and previously produced SYN2 families. Fifty-six of the 127 SYN2 families were created by using single plants from the F2 families included in Set 1. We will refer to these families as SYN2a and to the remaining 71 families as SYN2b. The production of SYN2 families followed the same procedure described for Set 1, which involved pooling and multiplication of the seed in isolated plots. Single plants were selected from the best F2 families by visual merits and by synchronous heading time.

Seeds for all families were produced at DLF Research Division in Store Heddinge (Denmark) and then shipped and sown in trials at other locations across Europe.

Sequence data were produced by genotyping-by-sequencing (GBS) (Elshire et al., 2011), which can be used in breeding populations to estimate genome-wide allele frequency profiles (Byrne et al., 2013). The GBS sampling and library preparation were performed according to Byrne et al. (2013) and Elshire et al. (2011). The plant material for the GBS data was obtained from remnant F2 seeds. Samples were digested using ApeKI, a methylation-sensitive restriction enzyme (to target the low copy fraction of the genome) and ligated to adaptors for sequencing and barcodes for individual family identification. This makes it possible to have separate data sets for each family and to estimate allele frequencies at each single-nucleotide polymorphism (SNP) location within each F2–SYN2 family. We prepared 32 libraries, each with up to 64 families, and sequenced each library on multiple lanes on an Illumina HiSeq2000 (single-end). The average number of reads obtained for each family after basic data filtering was ∼20 million. Data for each population were aligned against a draft sequence assembly (Byrne et al., 2015) and ∼1.8 million SNPs were identified with per-family sequencing depth at a SNP ranging from 1 to 250 (upper limit) reads per family. Few SNP positions had more than 60 reads and we suspect that these reads may be originating from plastid genomes or from highly repetitive regions not captured in the draft assembly. We therefore decided to discard all SNPs with average depth >60. Further, SNPs with allele frequencies <0.02 or >0.98 were removed. After this filtering on sequencing depth and frequency, 1,447,122 SNPs were available for analysis.

Phenotypic data were collected on family means during the standard breeding procedures of DLF A/S. Analyses focuses on various traits with different heritability:

  1. Crown rust resistance (CRR) measured by visual scoring during the period of maximum infection. The scale ranged from 1 (plant completely covered by rust) to 9 (no rust symptoms).

  2. Seed yield (SY), expressed in g m−2.

  3. 1000-kernel weight (TKW), which is the weight of 1000 seeds (g). Five technical replicates for each family were sampled, weighted, and then counted by digital imaging.

  4. Neutral detergent fiber (NDF), which is an indicator of the total fiber content (lignin, hemicellulose, and cellulose). It was measured in an ANKOM 2000 fiber analyzer according to the manufacturer’s instructions (Ankom Technology, 2013) and expressed as percentage of total dry matter content.

  5. Fructan (FR), which is the predominant type of storage carbohydrate in perennial ryegrass. It was measured as described in Bertram et al. (2010) and expressed as percentage of the total dry matter.

The CRR was recorded over 13 yr (from 2001 to 2013) across five locations (Table 1). All fields were organized in trials, which were further divided in subtrials and plots. In Les Alleuds (France) and Flevo Polder (Netherlands), families were sown during spring in small plots (0.5 by 4 m) and scored in late summer to autumn after every crown rust attack (1–3 times per year). Plots were cut between the different scoring time points to make them independent. In the other locations, the trait was scored on standard sward plots, in which families were sown during summer and farmed over two cropping seasons according to local management schemes (Table 1). A detailed description of the field design is given in Fè et al. (2015b). The other traits were measured in Store Heddinge (Denmark) between 2010 and 2013. Families were sown during spring in 1.5 by 10 m plots and farmed until seed harvesting (first half of August in the following year). The NDF and FR were measured on samples that were cut on the day of heading and allowed to dry in the field for 24 h. However, because of the high number of plots, some families (especially intermediate ones) were cut several days after heading date. Families were always sown in replicates, except for CRR in 2001, in Store Heddinge and Didbrook. Replicates were always sown in different trials. The experimental design was rather unbalanced. The level of unbalance was lower in the case of CRR, for which almost all families were sown at least in two different locations and in two different years. Data for TKW were measured only on one of the two replicates. The position of families within trial was randomized. However, randomization across trials was not always optimal, especially in the oldest experiments, where the field design was related to heading date, resulting in a certain degree of assortment between trials and the pedigree. A summary of the phenotypic data is presented in Table 2. For all traits, the number of phenotyped families and the number of locations, environments(location × year), and scores(location × year × scores) were recorded, along with the descriptive statistics (mean, standard deviation [SD], minimum, and maximum).

View Full Table | Close Full ViewTable 1.

Geographic positions, soil properties, and management conditions for seven experimental locations.

Location Position Soil Management
Flevo polder S Netherlands Loam Small plots
Les Alleuds W France Loam Small plots
Rennes NW France Light clay Conservation management (4–5 cuts per year)
Didbrook S England Clay and stony First year: conservation management (5 cuts); Second year: simulated grazing (7–9 cuts)
Store Heddinge SE Denmark Clay Conservation management (4–5 cuts per year); plus seed yield plots

View Full Table | Close Full ViewTable 2.

For all traits or set, number of phenotyped families, number of locations, environments (location × year interactions) and scores, mean, standard deviation (SD), minimum and maximum values.

Trait and set No. families No. locations No. location × year interactions No. scores Mean SD Min. value Max. value
Crown rust resistance
 F2 families 1790 5 24 49 5.35 1.85 1 9
 Synthetic families 66 5 10 22 5.64 1.53 1 9
Seed yield
 F2 families 1646 1 3 138.66 34.26 16.80 288.47
 Synthetic families 80 1 3 146.64 30.38 71.53 216.87
1000-kernel weight
 F2 families 1028 1 2 2.13 0.34 0.90 3.25
 Synthetic families 0
Neutral detergent fiber
 F2 families 1102 1 2 53.55 4.25 39.77 63.77
 Synthetic families 16 1 1 51.12 2.48 46.69 55.29
 F2 families 1544 1 3 7.63 4.07 0.35 29.42
 Synthetic families 32 1 2 4.87 2.90 1.93 14.95

Statistical Models and Estimation of Genetic Parameters

Data were analyzed by using different linear mixed models. Models were tested and compared using F-tests for the fixed part and Akaike tests for the random part. Genomic information was incorporated by genomic best linear unbiased prediction (GBLUP) (Habier et al., 2007; VanRaden, 2008). Each trait was analyzed by two different models: Model GE with the genomic effect as the only genetic effect and Model GP, which includes both PP effect and the genomic effects, allowing for split of genetic effect in two components: one among PPs and one within PPs and among F2 families. The models that had the best fit to the data are presented below.where y is the vector of observations; X is the design matrix of fixed factor; t is the vector of trials within locations and years (YLT); s is used only for CRR and refers to multiple scores performed within the same year; ts is the vector of trials within years, locations, and scores (YLST); d represents the difference between the date of heading and the date of harvest, expressed in number of days and divided in five classes (0–4, 5–8, 9–12, 13–16, > 16); td is the vector of d nested within YLTs; Zi are design matrices of random factors; i is a vector of breeding values ∼ N(0, Gσ2i), where G is the genomic relationship matrix; il is a vector of genotype × location interactions (accounts for the presence of replicates in certain fields) ∼ N(0, Iσ2il); ily is a vector of genotype × location × year interactions ∼ N(0, Iσ2ily); c is a vector of spatial effects within YLTs ∼ N(0, Iσ2c); p is a vector of the originating PPs ∼ N(0, Pσ2p), where P is a relationship matrix between PPs; ply is a vector of originating PPs nested within locations and years ∼ N(0, Iσ2ply); pp is the vector of interaction between PPs ∼ N(0, Iσ2pp); e is a vector of random residuals ∼ N(0, Iσ2e). In the models, CRR effects accounting for the interaction between single scores and the other factors are also present: ilys is a vector of genotype × location × year × score interactions ∼ N(0, Iσ2ilys); plys is the vector of interactions between PPs, location, year, and scores ∼ N(0, Iσ2plys); cs is the vector of spatial effects within YLSTs ∼ N(0, Iσ2cs). Finally, o is the vector of plots within YLTs across different scores ∼ N(0, Iσ2o). An additional effect of breeding values ∼ N(0, Iσ2i) was included to check the presence of nonadditive effect not captured by the genomic effect. This effect was never significantly different from zero and was therefore discarded from the analyses. The decision to use t as fixed effect was made because of the partial confounding between trials and the pedigree to avoid overestimation of the genetic variance (Henderson, 1984). Such a confounding was due to the fact that families were often sown according to earliness or according to the year in which the breeding material was originated. Design matrices for p, ply, and plys were built by accounting for the presence of multiple PPs as shown in Supplemental Figure S1: in each row, the numbers indicate the expected probability for each allele to come from a certain PP. The numbers on each row sum up to two, as each locus has two alleles. The G and P matrices were based on all SNP markers (after filtering for frequency and sequencing depth) and computed according to the first method described by VanRaden (2008). As genotypes were expressed as allele frequencies, and as they could only describe the genomic variance among families, the formula was modified according to Ashraf et al. (2014). Allelic frequencies for PPs were estimated as in Fè et al. (2015a). Off-diagonals of G (both within and between sets) and P, respectively, showing the relationship among families and PPs are reported in Supplemental Figure S2.

Variance components were estimated by restricted maximum likelihood method using the software DMU (Jensen et al., 1997; Madsen and Jensen, 2013) and can be interpreted as: σ2i is the additive genetic variance among families which is explained by the genomic effect tested across locations and years; σ2il is the genotype × location variance assumed to be constant over locations; σ2ily is the genotype × location × year variance; σ2ilys is the G × E variance among F2 families from repeated scores (nested within locations and years); σ2c is the variance among subtrials within YLTs; σ2cs is the variance among subtrials within YLSTs; σ2o is the environmental variance within plots and across scores; σ2p is the variance among PPs across locations and years; σ2ply is the G × E variance from PPs; σ2plys is the G × E variance from repeated scoring among PPs within locations and years; σ2pp is the variance among PPs combinations; and σ2e is the variance of residuals, which includes measurement errors and microenvironment effects within plots.

Genetic population parameters were calculated based on the additive biallelic infinitesimal model proposed by Ashraf et al. (2014) and already used in Fè et al. (2015a; 2015b). The original model ignores the relationship between PPs, which would cause inbreeding between the F1 families (Ashraf et al., 2014) and, as a consequence of that, changes in frequencies (and in variances) both among and within PPs and F2 families. That was taken into account by the use of P-matrix, for the among-PPs component and by the G-matrix, for the between PPs and among F2 families component. Assumptions are the same as stated in Fè et al. (2015b) except for the relatedness among PPs. Narrow-sense heritabilities (h2n; across environments, scores, and PPs) were calculated for both Model GE and Model GP. Here, only equations for SY are displayed, as formulas for the other traits can be easily written by adding or removing variance components from the denominator:where and are the average diagonals of G and P, respectively. Both are defined as the heritability of the family mean measured on a single plot. The component σ2p was added twice, as each row in the respective design matrix sums up to two (Supplemental Figure S1).

Cross-Validation Schemes

Prediction accuracies were first computed among the entire set of F2 families using GBLUP. Phenotypes were corrected for the fixed effect by running the model on the whole dataset. The GEBVs were then estimated by cross-validation by deleting phenotypes according to two different schemes that test for different hypothesis: (i) k-fold (k = 100), which tests predictions in case of presence of related individuals in the training and in the validation set, leaving out 1% of the families in each round, in random order; and (ii) pp-fold, which tests predictions in case of absence of related individuals in the training and in the validation set, estimating all the families originated by a certain PPs combination, after having left out everything that had at least a PP in common. As pp-fold implied a greater reduction in terms of training population than k-fold, a pp-like strategy was also tested to ensure a proper comparison. The pp-like strategy exactly replicated the cross-validation scheme used in pp-fold leaving out the same number of families in each cycle but choosing them at random instead than based on the PPs. This strategy ensures that the same size of training population is used in both schemes. All strategies were tested for both Model GE and Model GP. As the number of phenotyped families was different among traits and sets, to allow a proper comparison between accuracies, all models and strategies were also run on different sets with reduced population size consisting of randomly chosen F2 families. Such analyses were repeated 100 times, each time using a different set of randomly chosen F2 families, and the average predictive abilities and bias were calculated. The average accuracies (defined below) are reported in the result section. Finally, all F2 families were used to predict phenotypes for the three sets of synthetic families (SYN2, SYN2a, and SYN2b).

Predictive abilities and accuracies were calculated as described in Fè et al. (2015a). Briefly, accuracy, defined as the correlation between true breeding values and GEBVs (ρg,ĝ), was estimated with the following equation (Su et al., 2012):where ӯf represents the mean phenotypes corrected for the fixed effect and the numerator is the correlation between ӯf and GEBVs (predictive ability). The denominator represents the expected correlation between breeding values and corrected phenotypes and can be calculated from the following equation (Crossa et al., 2010):where σ2g is the genomic variance, and n is the number of replicates per each genotype. That is equivalent to the square root of the heritability based on several observations and represents the upper limit for accuracy of predicting phenotypes. However, this equation refers to a simplified model having only genomic and residual variances. In this paper, the formula needs to account for the other random factors. Furthermore, as traits were analyzed with different models, different formulas were also used to compute ρӯf,g. Here, only equations for SY are reported, as formulas for other traits can be easily derived:where is the average number of replicates per each family across all fields (nplots/nfamilies) and is the average number of environments per family (nily/nfamilies). Bias of the predictions was investigated by regressing ӯf on the breeding value estimates:

A regression coefficient (b) of 1 indicates no bias in the estimation of GEBVs. Comparison between Model GE and Model GP was assessed by correlating the respective solutions for GEBVs, and differences were tested using Hotelling–Williams test (Dunn and Clark, 1971) on ρӯf,ĝ using the R script developed by Christensen et al. (2012). In case the models gave a different correction for the fixed effect, ӯf was taken from the model that had the best fit to the data. Presence of reranking was tested by ranking families based on GEBVs from Model GE and Model GP, selecting the best 10%, and counting the number of families in common between the two different models.


Variance Components and Heritabilities

For each trait, the total phenotypic variance for all F2 families is shown in Table 3 together with the heritabilities and their standard errors (SE). Heritabilities from Model GE and GP were not significantly different. For this reason, only results from Model GP are shown. The percentage of each variance component over the total phenotypic variance is displayed in Fig. 1 (numerical values are available in the Supplemental Material). Narrow-sense heritabilities ranged between 0.40 and 0.55 for SY, TKW, NDF, and FR. Heritability for CRR (0.27) was lower than for the other traits investigated.

View Full Table | Close Full ViewTable 3.

Phenotypic variance (σ2P) of F2 families and narrow-sense heritability (h2) across environment and scores (with SE) for all traits.

Trait σ2P h2n
Crown rust resistance 1.914 (0.18) 0.27 (0.02)
Seed yield 477.2 (58.3) 0.42 (0.05)
1000-kernel weight 0.053 (0.01) 0.47 (0.14)
Neutral detergent fiber 10.68 (1.24) 0.54 (0.07)
Fructan 5.648 (0.58) 0.40 (0.06)
Fig. 1.
Fig. 1.

Percentage of variance components over the total phenotypic variance for all traits. CRR, crown rust resistance; SY, seed yield; TKW, 1000-kernel weight; NDF, neutral detergent fiber; FR, fructan. See text for variance component definitions. i, additive genetic variance within parent populations (PPs); p, additive genetic variance among PPs; il, genome × location variance; ily, genome × location × year variance; ply, genome × location × year variance among PPs; ilys, genome × location × year × score variance; plys, genome × location × year × score variance among PPs; pp, variance among PPs combinations; c, variance among sub-trials within YLTs; cs, variance among sub-trials within YLST; o, environmental variance within plots and across scores; e, residual variance.


In Fig. 1 and Supplemental Table S1 it is possible to distinguish between the amount of additive genomic variance that was present among PPs (σ2p) and the one within PPs (σ2g). Variance from PP (σ2p) was estimated to be 39 and 29% of the genomic variance in NDF and TKW, respectively. It was lower in SY (18%) and FR (16%) and very low in CRR, accounting for only 9% of the total genetic variance and the 2.5% of the phenotypic variance. The G × E effects turned out to be low in SY, representing only 8.6% of the phenotypic variance, and large in CRR, where they accounted for 30% of the total variation. For CRR, the strongest effect (half of the total G × E variation) was from the genome (+PPs) × location × year × scores component. The other half of G × E was divided between genome × location (18%) and genome (+PPs) × location × year (32%).

Relationship Between Sets and Cross-Validations Schemes

Relationship between sets can be inferred from the boxplots of genomic relationships reported in Supplemental Fig. S2. Genomic relationships between F2 families ranged from −0.2 to almost 1.0, indicating the presence of high relationship within this set. Off diagonals between F2 and SYN2 families showed a similar distribution but with less extreme values. As expected, relationships between F2 and SYN2a families were greater than the ones between F2 and SYN2b but not dramatically. Results from cross-validation are displayed in Table 4 and in Fig. 2. Probably as a result of the small number of families, accuracies for the three sets of SYN2 families (SYN2, SYN2a, and SYN2b) were not statistically different from each other. Therefore, only results for SYN2 are reported. All predictive abilities (and accuracies) were significantly different from zero (P < 0.001). Scheme pp-like, depending on the model and trait, was either equal to k-fold or gave slightly lower accuracies with the highest drop being 4%. Therefore, results for that scheme were not reported. For Model GE, when related material was present in the training population (k-fold scheme), all accuracies were >0.50. The best estimates were registered for SY in both the full F2 sets (Table 4) and in the sets with reduced population size (Fig. 2). The regression of corrected phenotypes on GEBVs did not indicate the presence of significant bias for TKW and FR. For the other traits, regression coefficients were between 1.12 and 1.18, indicating an upward bias in the variance of the genomic predictions. As expected, for all traits, increasing the number of individuals used to train the model led to an increase in accuracy (Fig. 2). This effect was mostly pronounced in population sizes up to 500 to 750 families, depending on the trait, except for FR, for which accuracies kept increasing significantly also with larger population sizes. Predictions for population sizes of <175 families are not shown because they were affected by very large SE and not indicative of any trend.

View Full Table | Close Full ViewTable 4.

For all traits and models, predictive ability (ρӯf;ĝ), accuracy (ρg;ĝ), and bias (with SE), for the k-fold and pp-fold cross-validation, and for the prediction of synthetic (SYN) families.

Trait and model† k-fold
F2 → SYN2
ρӯf;ĝ ρg;ĝ Bias (SE) ρӯf;ĝ ρg;ĝ Bias (SE) ρӯf;ĝ ρg;ĝ Bias (SE)
Crown rust resistance
 GE 0.58 0.68 1.18 (0.04) 0.45 0.53 1.38 (0.07) 0.54 0.68 0.90 (0.17)
 GP*** 0.57 0.69 1.27 (0.04) 0.46*** 0.56 1.47 (0.07) 0.48 0.61 0.96 (0.21)
Seed yield
 GE 0.56 0.73 1.12 (0.04) 0.48 0.63 1.13 (0.05) 0.46 0.50 0.87 (0.22)
 GP*** 0.58** 0.75 1.09 (0.04) 0.48 0.62 1.17 (0.05) 0.45 0.49 0.83 (0.19)
1000-kernal weight
 GE 0.37 0.51 1.06 (0.08) 0.31 0.42 1.16 (0.11)
 GP** 0.44* 0.65 0.99 (0.06) 0.37*** 0.54 1.68 (0.13)
Neutral detergent fiber
 GE 0.56 0.69 1.13 (0.05) 0.48 0.59 1.33 (0.07)
 GP*** 0.68*** 0.81 1.02 (0.03) 0.56* 0.68 1.70 (0.08)
 GE 0.37 0.50 1.05 (0.07) 0.27 0.36 0.99 (0.09) 0.48 0.53 1.62 (0.54)
 GP*** 0.45*** 0.60 0.99 (0.05) 0.31*** 0.41 1.13 (0.09) 0.47 0.51 1.37 (0.48)
*Significant at the 0.05 probability level.
**Significant at the 0.01 probability level.
***Significant at the 0.001 probability level.
Asterisks indicate significant of differences from Akaike test between Model GE and GP.
Asterisks for ρӯf;ĝ indicate the significance of differences between ρӯf;ĝ from Model GE and GP from Hotelling–Williams test.
Fig. 2.
Fig. 2.

Accuracies in the reduced set of F2 families for all traits in the k-fold cross-validation strategy for Model GE (gray line) and Model GP (black line). CRR, crown rust resistance; SY, seed yield; TKW, 1000-kernel weight; NDF, neutral detergent fiber; FR, fructan.


Accuracies, when predicting from unrelated families (pp-fold scheme), were significantly worse than the ones obtained for k-fold and pp-like (8–16% lower). Bias was generally higher but not significantly different from what was found for the k-fold scheme.

Finally, prediction of multiparental SYN2 families from biparental F2 families (Table 4) showed high accuracies for SY (0.50), FR (0.53), and CRR (0.68). Because of the low number of families (only 16) it was not possible to get significant accuracies for SYN2 families in NDF. Bias for SYN2 families’ prediction was not significant.

Implementation of Parental Populations Effect

The inclusion of PPs effect generally improved predictions. The Akaike test showed Model GP to have a best fit to the data in all traits. When breeding values were predicted, also from related families (k-fold cross-validation), predictive ability (and accuracy) was significantly higher for all traits, except for CRR. However, for CRR, as shown in Fig. 2, Model GP gave better predictions in case of smaller training population sizes. The increase in predictive abilities was moderate for SY (+0.02) and very high for NDF, FR, and TKW (∼ +10). However, correlating GEBVs with the same ӯf showed the actual increase in correlation to be ≤0.04. For TKW, the p-value was higher because of the smaller amount of data and the larger error. After PPs implementation, the best accuracy was obtained for NDF for both large and small training population size. Including PPs also decreased the bias in all traits (especially for NDF) with the exception of CRR.

When no related families were present in training and validation set (pp-fold scheme), the implementation of PPs had smaller effect on accuracies for NDF and no effect for SY. On the other hand, it caused a significant improvement in the case of CRR, TKW, and FR. No positive effect was found for bias; regression coefficients (b) were significantly >1 for all traits.

Furthermore, including PPs in the prediction models generally showed significant effects on the ranking of families based on GEBVs (Table 5). Correlations between GEBVs from the two models (GE and GP) were always significantly different from 1 in the k-fold scheme (0.93–0.97). Values were higher in the pp-fold scheme (0.96–1.0), and not significantly different from 1 in the case of FR and SY. In these cases, the percentage of common families within the first 10% of the ranks (genotypes that will be selected for further steps of the breeding program) was around 90%. For correlations ∼0.93, the percentage was between 65 and 75%.

View Full Table | Close Full ViewTable 5.

Correlations between genomic estimated breeding values (GEBVs) from Models GE and GP and percentage of common F2 families in the best 10% for k-fold and pp-fold cross-validation.

Trait k-fold
ρGEBV(GE);GEBV(GP) Common F2 ρGEBV(GE);GEBV(GP) Common F2
Crown rust resistance 0.97 0.84 0.96 0.81
Seed yield 0.97 0.82 0.99 0.91
1000-kernel weight 0.93 0.76 0.97 0.82
Neutral detergent fiber 0.94 0.72 0.97 0.83
Fructan 0.94 0.66 0.99 0.89

Finally, PPs implementation did not result in any significant improvement for prediction of the SYN2 families in terms of predictive ability and accuracy nor in terms of bias.


The present study demonstrates the potential to implement GS in breeding programs for perennial ryegrass, confirming expectations from previous theoretical studies (Hayes et al., 2013).

Variance Components and Heritabilities

A considerable amount of additive genetic variance (σ2g + σ2p) was found in all traits. A large contribution to the total variance for CRR also came from G × E interactions, especially from the genotype × environment × scoring component (σ2ilys + σ2plys). This observation agrees with the current state of knowledge on this trait, which is believed to be determined by both quantitative genes, mainly involved in slow-rusting resistance, and major genes, likely responsible for race-specific resistance (Dracatos et al., 2010). The size of σ2ilys relative to σ2ily and σ2il indicates that the level of G × E within fields and years was greater than the level of variation between locations and environments. This may not be surprising, as CRR was shown to be highly affected by a wide range of environmental and physiological factors such as temperature, light, relative humidity, leaf surface topography, and health and susceptibility of the host plant (Dracatos et al., 2010). The pathogen also seems to display a significant genetic variation and constantly migrates and rapidly evolves through mutations, recombination, and crosses between different races (Dracatos et al., 2010).

Regarding FR, concerns may be raised by the significance of the factor d (difference between heading date and actual harvest date), which also showed to vary across trials. This fact may be related with the high variability over time of FR content, which was also shown to be subject to diurnal regulation (Longland et al., 1999). We did not observe similar effects on NDF, which is another forage quality trait. Neutral detergent fiber accounts for most the structural parts of the plant (cellulose, hemicellulose, and lignin), and previous investigations revealed that while NDF generally increase over time, this effect is most pronounced for early heading varieties (data not shown).

The level σ2p compared with σ2g showed the total genetic variance to be greater within PPs than among PPs. This observation was previously reported by Fè et al. (2015b) on a subset of this plant material and was true for a wide set of traits with the exception of heading date. This can be explained by the fact that families and plants are usually crossed with genetic material that has similar heading date. This assortative mating was never performed for NDF, FR, SY, or TKW, so the degree of variance between PPs may be due to a certain degree of correlation with heading date. For NDF and FR, this correlation has been documented in previous studies (Humphreys, 1991; Sampoux et al., 2011).

Cross-Validations: General Trends

Predictive abilities and accuracies were high for all traits. The GEBV estimations were better when families were predicted from related breeding material (scheme k-fold and pp-like) but were significantly worse when related individuals were left out from the training population (scheme pp-fold). In a standard breeding program, where new PPs are regularly added to an existing genetic pool, the situation is likely to be intermediate between k-fold and pp-fold. Therefore, to improve predictions, it may be a good idea to always phenotype the newly introduced PPs or a part of their offspring.

Analysis on different population sizes showed that for almost all traits under the present breeding scheme, at least 500 families are needed in the training population. The same outcome was previously reported on the same breeding material for heading date (Fè et al., 2015a). Therefore, it may be that the relation between training population size and accuracy depends mainly on the population structure rather than on the nature of the trait. The exception of FR may be related to the higher number of fixed factors in the models shown in Eq. [9] and [10], which limits the statistical power of the predictions.

Predictions of multiparental SYN2 families were possible only for CRR, SY, FR, and for a small set of families because of lack of records. Estimates for the two traits cannot directly be compared with each other because the phenotypic dataset for CRR and SY only had a few families in common. Accuracies for SYN2 families were remarkably high and showed good possibilities for predicting breeding values. That was true even for SYN2 families that originated by breeding material not included in Set 1. This fact is likely to be related to a certain degree of relationship between F2 and SYN2b families (Supplemental Fig. S2), which probably arises from relationship between PPs.

The bias found among F2 families indicated a general tendency in predicting GEBVs with a too small variance, which may be related to the partial confounding between the trials and the pedigree. Such unbalance can cause a small part of the genetic effect to be captured by the fixed part of the model.

Implementation of Parental Populations Effect

The inclusion of the PPs effect always improved the fit to the data, but its role in improving predictions was less clear. The effect was of course higher in traits characterized by a large variance among PPs (σ2p), which was the case for NDF and TKW. The increase in accuracy mainly depended on two factors: (i) actual improvement in breeding value estimation (quantifiable in an increase in predictive ability up to 4%), and (ii) better fit to the data and, consequently, better correction for the fixed effect. The second point is especially important for NDF, and it may be related to the field design in which PPs and trial within locations and years effects were confounded to some degree. The lack of improvements when implementing PPs in the prediction of SYN2 families may be explained by the presence of a high number of PPs used to generate the SYN2 families.

Regarding bias, the reduction brought by Model GP in the k-fold scheme, compared with Model GE, may be explained by a better separation between the genetic and the spatial effect from the introduction of the effect of PPs (p). Furthermore, in many cases, including PPs in the model resulted in significant changes in the ranking of families based on GEBVs and, because of that, also in the list of families that will be selected for the following steps of the breeding program.

Cross-Validations: Trends for Each Trait

Results reported for SY were very encouraging, showing prediction accuracies of 0.75 when predicting among all breeding material and between 0.50 and 0.60 when predicting from families without PPs in common. The outcome is very positive also when compared with accuracies reported for other crops. So far, the trait has been widely investigated only in wheat, where accuracies where found between 0.2 and 0.79, and maize, which showed accuracies ranging from 0.42 to 0.5 (Lin et al., 2014). The high potential in breeding value prediction, together with a fairly high heritability, makes it a very promising target trait for GS implementation. A significant improvement for this trait will be extremely beneficial for breeders, as a poor seed production represents one of the most limiting factors in several breeding schemes. However, when improving SY, attention must be paid on its correlations with forage yield and digestibility. These correlations have been shown to be somewhat negative but possible to break in an efficient breeding scheme (Studer et al., 2008), especially if focusing on SY components such as tiller numbers or proportion of reproductive tillers.

Regarding CRR, to our knowledge, this paper represents the first example of prediction based on genomic information. However, articles were recently published on wheat focusing on similar diseases (yellow rust and stem rust) and showing overall good perspectives (Ornella et al., 2012; Crossa et al., 2014). Predictive abilities were variable ranging from zero to 0.80 depending on the environment and the nature of the breeding material. Predictive ability across populations was 0.62 (Crossa et al., 2014), which is similar to what was found in this paper for the k-fold scheme. In this study, accuracies were relatively high even in the pp-fold scheme and in predicting SYN2 families from F2 families despite the relatively low heritability. Overall, this study showed very good potential for predicting and improving the genomic effect acting across environments and scores; however, the issues concerning G × E still need further attention. Because of the small σ2il and the large σ2ilys, rather than focusing on breeding varieties for specific environment, it may be fruitful to test models that also include weather and physiological conditions to have a better model for G × E interactions (e.g., Technow et al., 2015) and also develop models that account for different pathogen races.

Positive results were also achieved for NDF and, to a lesser extent, for TKW and FR. In this case, comparison with other species is complicated, as the number of publications is limited. Thousand-kernel weight was investigated in wheat giving accuracies between ∼0.30 (Poland et al., 2012) and ∼0.80 (Lado et al., 2013). Regarding fiber content, a study on sugarcane (Saccharum officinarum L.) (Gouy et al., 2013) reported predictive abilities up to 0.40 for acid detergent fiber, which is one component of NDF (Goering and Van Soest, 1970). For these traits, phenotyping the same families in another location may be useful to better understand the role of G × E interactions. Furthermore, for FR, in further studies, it will be very important to sample all families, ensuring similar sampling conditions (e.g., developmental stage and harvest hour).


This study demonstrates that the implementation of GS in grass breeding is possible and presents an opportunity to make significant gains for economically important traits. We found a significant level of genetic variance for all traits studied, a very large proportion of which could be traced by genomic information from genotyping assay.

Of concern is the presence of significant G × E interaction for some traits especially in resistance to crown rust. On the other hand, accuracies were high both among F2 families and between biparental F2 families and multiparental synthetics, demonstrating a potential for implementing genomic prediction at different stages of the breeding cycle and using different types of breeding material. Further studies are needed on determining the stages in which such implementation would be most effective and whether it would be advantageous to restructure the current breeding scheme.

Accuracies were significantly higher when training and validation set contained related individuals. The study also indicates an advantage in implementing the effect of PPs in the models for prediction compared with a model with genomic effects of F2 alone. However, such advantage seems to depend on the mating structure, as it was shown to be low when families were originated from many PPs.


The project was financed by the Danish Ministry of Education through the Council for Industrial Ph.D. Education (11-109967), by the Danish Ministry of Food, Agriculture, and Fisheries through The Law of Innovation (3412-09-02602) and the Green Development and Demonstration Program (GUDP; 3405-11-0241), and by DLF A/S.

We also acknowledge Adrian Czaban and Stephan Hentrup, who collaborated on the preparation of GBS libraries, working at Aarhus University, Department of Molecular Biology and Genetics, Crop Genetics and Biotechnology.





Be the first to comment.

Please log in to post a comment.
*Society members, certified professionals, and authors are permitted to comment.