Divergent Genomic Selection for Herbage Accumulation and Days-To-Heading in Perennial Ryegrass

: Increasing the rate of genetic gain for dry matter (DM) yield in perennial ryegrass ( Lolium perenne L.), which is a key source of nutrition for ruminants in temperate environments, is an important goal for breeders. Genomic selection (GS) is a strategy used to improve genetic gain by using molecular marker information to predict breeding values in selection candidates. An empirical assessment of GS for herbage accumulation (HA; proxy for DM yield) and days-to-heading (DTH) was completed by using existing genomic prediction models to conduct one cycle of divergent GS in four selection populations (Pop I G1 and G3; Pop III G1 and G3), for each trait. G1 populations were the offspring of the training set and G3 populations were two generations further on from that. The HA of the High GEBV selection group (SG) progenies, averaged across all four populations, was 28% higher ( p < 0.05) than Low GEBV SGs when assessed in the target environment, while it did not differ significantly in a second environment. Divergence was greater in Pop I (43%–65%) than Pop III (10%–16%) and the selection response was higher in G1 than in G3. Divergent GS for DTH also produced significant ( p < 0.05) differences between High and Low GEBV SGs in G1 populations (+6.3 to 9.1 days; 31%–61%) and smaller, non-significant ( p > 0.05) responses in G3. This study shows that genomic prediction models, trained from a small, composite reference set, can be used to improve traits with contrasting genetic architectures in perennial ryegrass. The results highlight the importance of target environment selection for training models, as well as the influence of relatedness between the training set and selection populations. from synthetic populations assessed as plots in a commercial breeding program These studies have demonstrated the potential of GS as a tool to accelerate improvement for agronomic traits in ryegrass, as well as indicating ways forward for the implementation of GS in breeding programs.


Introduction
Pastures based on perennial ryegrass (Lolium perenne L.) are the principal source of nutrition for ruminant livestock in temperate environments across the globe. The harvestable dry matter yield (DM yield) is regarded as an important factor influencing livestock production per hectare [1]. The profitability of dairy farms in countries such as New Zealand and Australia, and the international competitiveness of these industries, are closely linked to the consumption of homegrown forage [2] because it is a low-cost feed source. The contribution of DM yield to farm profitability is formally The relative infancy of GS in forage plant breeding, and the time required to evaluate long-term traits such as DM yield, mean that there are currently few reports on the impact of GS on realized performance. The purpose of this study was to determine if the application of divergent genomic selection using genomic prediction models trained for herbage accumulation (HA; [20], a proxy for DM yield potential), would elicit a selection response from breeding populations with differing levels of relatedness to the original training set. As a comparison, independent genomic selections were also made for days-to-heading (DTH), which is a genetically simpler trait of higher heritability, using a prediction model developed with the same training set.

Genotyping and Determination of GEBVs
Models for the genomic prediction of herbage accumulation (HA) and days-to-heading (DTH) were developed previously, using data from a training set composed of five discrete breeding populations (Pop I to V) [20]. In that study, prediction models were trained using genotyping-bysequencing (GBS) SNP genotypes from the maternal parent plants in the training set, with phenotypic data acquired from HS families of the respective maternal parents by evaluating sown plots under field conditions at multiple locations. In the current study, the HA and DTH prediction models developed were applied for divergent genomic selection in each of four synthetic selection populations (Table 1) sourced from the Grasslands Innovation Ltd breeding program. In all cases, the selection populations were descended from two of the corresponding training set populations (Pop I or Pop III) that had been used to derive the HA and DTH prediction models ( Figure S1). The selection populations were either the direct offspring (Pop I G1 and Pop III G1) of the training set generation or were from a more advanced generation (Pop I G3 and Pop III G3). The G3 selection populations had arisen because, while the training set trials described in [20] were being conducted, breeding from training set populations had continued separately in the Grasslands Innovation Ltd breeding program (for non-DM yield traits). The G3 selection populations are therefore advanced generations derived after two further cycles of conventional selection and recombination. GEBVs were derived for a total of 970 individuals from the four selection populations (n = 177-280 individuals per population; Table 1). For each of Pop I G1 and Pop III G1, there was a known HS family structure. For those populations, seeds for germination, genotyping, and subsequent GEBV determination were taken from six randomly-chosen HS families within each population. For Pop I G3 and Pop III G3, there was no recognised family structure and for those populations, seeds for germination, genotyping, and GEBV determination were chosen randomly from a bulk sample. Seeds for all but Pop III G3 were germinated and grown under standard greenhouse conditions for approximately four weeks, until two tillers had emerged, before sampling tissue for DNA analysis. For Pop III G3, pseudostem tissue was supplied directly by Grasslands Innovation Ltd and provided for genotyping. DNA was isolated from 50-100 mg of harvested pseudostem tissue using a semi-automated method for the extraction of sequencing-quality high molecular weight DNA from plants [27]. Each seedling was also subject to an immuno-detection blot [28], to confirm the presence of fungal endophyte (Epichloë festucae var lolii).
Methods for GBS library development, the sequencing of GBS libraries, sequence data processing, genotype calling, and genomic relationship matrix (GRM) estimates were as described in [20]. Briefly, 96-plex GBS libraries were generated for the 970 seedlings by digesting 100 ng of DNA with ApeKI (New England Biolabs, Ipswich, MA, USA) and ligating to a unique barcoded adapter and a common adapter (99 ng). Libraries were size-selected for fragments of 193-313 bp using a Pippin Prep TM DNA size selector (Sage Science, Beverly, MA, USA). Each library was sequenced on one lane of an Illumina HiSeq 2500 (Illumina, San Diego, CA, USA) at AgResearch Invermay, New Zealand. Raw reads from the selection population GBS libraries were merged with those previously generated for the training set [20] and SNP calling was conducted jointly using the TASSEL 5.0 GBS pipeline [29] with a ryegrass reference genome [20,30]. After filtering for missing SNPs (50% maximum missing per site), minor allele frequency (>0.05), and read depth (>2), the allelic counts of the remaining biallelic SNPs were exported to KGD software (https://github.com/AgResearch/KGD) [31] and further filtering based on Hardy−Weinberg disequilibrium (Hardy-Weinberg disequilibrium >−0.05, estimated as the observed frequency of the reference allele homozygote minus its expected value) was applied within the software, as recommended by [31]. A cut-off value of >−0.05 was applied to filter SNPs potentially from duplicated or repetitive regions of the genome. KGD is an approach for kinship estimation which accounts for the (low) read depth of the GBS data and does not require imputation. The KGD method outputs a GRM (may be referred to as a KGD matrix hereafter). GEBVs of the selection candidates were derived by Genomic Best Linear Unbiased Prediction (GBLUP) for the mean HA under each of the five environments described in [20], including the Waikato standard and severe summer defoliation, Manawatu standard and severe summer defoliation, and Canterbury standard summer defoliation. GEBVs for DTH were derived using a model that had been developed for a prediction of the mean DTH using data from two environments (Waikato and Canterbury, the latter in the South Island of New Zealand) [20].

Selection and Polycross Isolations
For HA, divergent selections were made in each selection population based on a simple additive index of HA GEBVs (models described in [20]), weighted towards growth in the Waikato environment, which is considered most challenging for ryegrass growth [32,33]: (1) (1) Both DTH GEBV and HA index GEBV were calculated for 970 individuals, but, due to either a nil endophyte status or seedling mortality, selections were confined to a total of 795 individuals (Table  1). Prior to selection, 15 individuals were randomly taken from each population to create a control (Non-selection group). Within the remainder of each selection population, the top 10% (High GEBV group) and bottom 10% (Low GEBV group) were identified for each of DTH and HA index, and 15 individuals were randomly selected from within each group for poly-crossing to enable progeny testing. The same Non-selection group was used for both HA and DTH in each selection population, so a total of 20 polycrosses were completed (four populations × two traits x High GEBV, Low GEBV; plus four populations × Non-selection groups) under isolation during spring 2016 at Palmerston North, New Zealand. Maternal HS family seed was harvested from individual plants within each polycross and balanced bulks made by combining equal seed quantities from all HSs to create a synthetic population (Syn1) for each of the 20 selection groups (hereafter referred to as SGs).

Evalaution Trials
Field trials for both DTH and HA were established at AgResearch Ruakura in Waikato (37.78° S, 175.32° E) and AgResearch Grasslands in Manawatu (40.21° S, 175.37° E), New Zealand. In both trials, soil fertility levels were adjusted to ensure that nutrients were not limiting plant growth. Nitrogen was applied (15-30 kg N/ha) at each defoliation. Superphosphate fertilizer (8.8 kg P/ha) was applied in late autumn each year. Plants in both HA and DTH trials were defoliated by sheep grazing whenever they reached the 2-3 leaf stage of development, except (a) when harvests were taken in the HA trials, in which case the plants were defoliated manually (see below), and (b) in the DTH trials from 8 September to 31 December, when the trials were closed to grazing to allow flowering to proceed.
For DTH, 25 random plants per SG were transplanted into spaced plant field trials in May 2017, at 40 cm spacing between plants in a row-column design with 45 repeated clonal checks (three check plants × 15 clones). DTH was measured in every plant as the number of days after 22 October until five fully-emerged flower heads were observed, in both 2017 and 2018. For HA index selections, the experimental SGs and two check cultivars were direct-drilled as 2 m rows (0.6 g seed per row) in a row-column design with n = 4 replicates. Eight discrete harvests (Table S1) were completed at each site by manual cutting, drying, and weighing of herbage to determine grams of DM per 2 m row, following the methodology described in [20]. HA was determined as the mean g DM per row across all harvests.
For DTH and HA trials, data were analysed by a linear mixed model, using the variance component analysis procedure residual maximum likelihood (REML) option in DeltaGen software v 0.03 (http://agrubuntu.cloudapp.net/PlantBreedingTool/) [34]. For DTH repeated clonal checks, year, location, and SG were treated as fixed effects, while rows and columns were treated as random. For HA, repeated checks, location, harvest date, and SG were treated as fixed effects, while replicates, rows, and columns were treated as random. The final SG values were generated as best linear unbiased estimates (BLUEs) [35].

GBS Data and GEBVs
Genotypes were called jointly amongst training set and selection population individuals and used to construct a genomic relationship matrix (Table S2) based on 777 k SNPs. Selection populations showed close genetic alignment to their training set progenitors (Figure 1), but the G3 populations were more distal than the G1 populations, most noticeably for Pop III. This was reflected in the mean relatedness coefficients estimated for selection populations and their training set counterparts (Table  1), with relatedness to the training set generation declining by 7% from training to G1 and by 26% from training to G3, in both Pop I and Pop III. Within-population genomic diversity was assessed based on the mean relatedness coefficient amongst individuals within the population. This was reduced in G1 and G3 generations for both Pop I and Pop III relative to the training generation, with a greater reduction in G3 ( Table 1). The decline in genetic diversity in the G3 populations is likely a result of selection imposed (for non-DM yield traits) between the training and G3 generations in both sets of material, as previously observed in ryegrass populations under selection [36]. The decrease in diversity between training and G1 generations was presumably affected by selection being focused on only six of the 112-121 HS families represented in the training generations of Pop I and III.
GEBVs for five HA traits (mean HA under Waikato standard summer defoliation, Waikato severe summer defoliation, Manawatu standard and severe summer defoliation, and Canterbury standard summer defoliation) and DTH were successfully generated for 970 individuals, using previously defined KGD-GBLUP prediction models for those traits. At this point, a total of 175 individuals were discarded from the selection process due to mortality or because they returned a negative endophyte immunoblot result (Table 1). Following GEBV-directed selection, poly-crossing, and seed harvest, progeny HS families were bulked into Low GEBV, Non-selected, and High GEBV selection groups (SGs) for each selection population x trait combination (n = 20 SGs in total; Non-selected SG's were shared for the two traits). The SGs were subsequently placed in field evaluation trials for the measurement of DTH or HA. Polycrosses yielded between 29.4 and 143.9 g seed (mean = 78.1 g), approximately equivalent to 15,000 and 72,000 seeds (mean = 39,000), assuming a thousand-seed weight of 2 g.

Genomic Selection for DTH
Days-to-heading (DTH) was investigated in this study as a test trait characterized by a moderatehigh heritability [25,37], in contrast to the low-moderate heritability trait-DM yield [38,39]. Days-toheading was measured in the SGs over two consecutive years at both Manawatu and Waikato sites and DTH BLUE values were estimated for each SG (Low GEBV, High GEBV, and Non-selected) from each of the four selection populations, across sites and years. Assuming successful divergent genomic selection, the expected selection response is a ranking of High GEBV > Non-selected > Low GEBV. A trend of High GEBV > Low GEBV was observed in all of the selection populations, but the magnitude of the response varied amongst the populations ( Figure 2). However, Non-selected SGs did not consistently rank between High GEBV and Low GEBV SGs in three of the populations and this asymmetry in the selection response is discussed further below. Statistically significant (p < 0.05) differences between High and Low GEBV SGs were detected in both G1 populations (6.3-9.1 days; 31%-61%), while smaller, non-significant (p > 0.05) differences were observed in their G3 counterparts (4.1 days in both; 13%-19%). While the latter is suggestive of a multi-generational genomic selection response, the apparent decline from G1 (direct offspring of the training set) to G3 (advanced generation) does indicate a decay in the predictive ability of the model. This decay was apparently associated with an increased genetic distance from the training set to G1 to G3, inferred from the mean relatedness (Table 1). The accuracy of genomic prediction models has often been observed to decrease with an increased genetic distance or number of meioses from the training set [40][41][42][43]. This is likely exacerbated by the use of GBLUP as the statistical method for conducting genomic prediction, in conjunction with the small training set size (n = 517) [20]. The GBLUP method principally exploits genomic relationships, which decay faster over generations than the more generationally-persistent linkage disequilibrium (LD) between SNPs and quantitative trait loci (QTL) [43,44]. The diminished response observed in the G3 populations may also have been affected by the level of available genetic variation in these populations, which was apparently lower than that of their G1 counterparts, based on estimates of within-population relatedness (Table 1). In certain selection populations, when comparing Low and High GEBV SGs to the performance of the Non-selected SG, the selection response was asymmetric. For example, in Pop I G1, genomic selection was most effective in increasing DTH, with no significant reduction in DTH in the Low GEBV relative to the Non-selected SG. The reverse was observed for Pop III G1, with no significant increase in DTH in the High GEBV relative to the Non-selected SG ( Figure 2). This may reflect differences between Pop I and Pop III in terms of heterozygosity or fixation at key QTL that affect late and early heading, respectively. This pattern of response was most pronounced when considering data from the Manawatu site alone, while the pattern was present, but diminished, at Waikato ( Figure S2). Waikato was one of the two locations in which the DTH genomic prediction model was trained, whereas Manawatu was not [20].
The lack of a significant selection response towards later heading in Pop III G1 is also interesting when placed alongside the fact that the mean DTH of the Pop III G3 population was significantly (p < 0.05) higher compared to its Pop III G1 counterpart (+6.8 days between Non-selected SGs; Figure 2). The higher DTH in Pop III G3 is consistent with the fact that phenotypic recurrent selection for later heading had been undertaken between G1 and G3 in the conventional breeding program (Dr Richard George, pers comm). This indicates that genetic variance for later heading must be present in Pop III G1, but for some reason, it was not captured by the genomic prediction model and therefore, there was no shift to a higher DTH in Pop III G1 when GS was applied. This may be a consequence of the relatively small, multi-population training set used to develop the prediction model-later-heading families from Pop III may not have been sufficiently represented in the training set. It may also be that large-effect QTL for the heading date were not in LD or linkage with SNPs in the training set and therefore not captured in the prediction model, although this seems less likely given the high density of SNPs used in the genomic prediction model. Although it has yet to be adequately demonstrated using real data [25], the use of non-GBLUP statistical prediction approaches, such as Bayesian models, which assume non-uniform SNP-specific variances, may be more effective for traits such as DTH, for which the genetic architecture is understood to include large-effect QTL [45]. The phenotypic selection response observed in Pop III G3 may also have been influenced by non-additive genetic effects [25] that were not accounted for in a genomic prediction model based on additive genetic variance.

Genomic Selection for HA
Herbage accumulation (HA) harvests were completed at both locations, at seven (Waikato) or eight (Manawatu) seasonal time points, between January and October 2019 (Table S1). The timing of harvests was not consistently aligned between sites and therefore, the assessment of the genomic selection response here focuses on HA BLUEs estimated for each SG (Low GEBV, High GEBV, and Non-selected) across all harvests, within the respective locations. This was completed for each selection population individually, as well as for the four populations combined (All Pop).
Based on the All Pop analysis, SGs performed better overall in the Manawatu environment compared with Waikato, based on the mean HA at the respective sites ( Figure 3, Table 2). The response to divergent GS for HA also showed an apparent effect of the environment. In the All Pop analysis, significant (p < 0.05) differences in HA amongst High GEBV, Low GEBV, and Non-selected SGs were detected at the Waikato site, but not Manawatu ( Figure 3, Table 2), and the expected ranking of High GEBV > Non-selected > Low GEBV was observed at Waikato only. As with DTH, there was an asymmetric nature evident in the genomic selection response at Waikato, albeit considerably influenced by one of the four selection populations-Pop I G1 (see below). At Waikato, the High GEBV SGs performed significantly (p < 0.05) better than the Low GEBV SGs (+28%), and the Non-selected SG was also significantly (p < 0.05) higher than the Low GEBV (+20%) ( Table 2). However, the difference between the High GEBV and Non-selected SGs was small and non-significant (+7%; p > 0.05). In addition to the potential underlying causes outlined for DTH above, asymmetry in the selection response for HA may also have been influenced by non-additive genetic effects, such as heterosis, which were not accounted for in the genomic prediction model, but which may be significant for the Syn1 generation studied here.
At Waikato, the expected response to selection was also observed when the four populations were considered individually ( Figure 3, Table 2). Strong responses occurred in both Pop I G1 and Pop I G3, with significant (p < 0.05) HA increases of 65% and 43% between the Low and High GEBV SGs of the respective populations. As noted above, the selection response in Pop I G1 was highly asymmetric, with no discernable difference between the High GEBV and Non-selected SGs. In Pop I G3, there was a similarly asymmetric selection response, but the Non-selected SG was more clearly intermediate compared to the Low and High GEBV groups, despite only being significantly different (p < 0.05) to the Low GEBV SG. Overall, these data indicate a stronger effect of genomic selection for reducing DM yield than for increasing DM yield. The expected ranking of SGs was also observed for both of the Pop III selection populations at Waikato, but differences amongst High GEBV, Low GEBV, and Non-selected SGs were not statistically significant (p > 0.05). At Manawatu, a directional response to GS was only observed in one of the four populations-Pop I G1 (Figure 3)-although differences amongst the SGs in that population were statistically non-significant (p > 0.05).  Our research is aligned to the paradigm of breeding for DM yield for broad adaptation across agricultural environments. However, a unified multi-environment HA genomic prediction model was not available for testing in this experiment, because at the training stage across-environment, adjusted means for HA were not successfully generated, largely due to the influence of significant genotype-by-environment interactions [20]. Instead, we conducted GS using a weighted index of GEBVs for HA performance in each of the five individual training environments (Equation (1)) and this index was substantially weighted (60%) towards Waikato environments. The difference in the selection response observed between the two locations, with a significant response seen in Waikato but not Manawatu, most likely reflects the weighting of the GEBV index and emphasizes the importance of phenotyping the training set within the target environment(s) [46]. In [20], it was noted that developing a model for the genomic prediction of HA across all New Zealand environments may be difficult to achieve, due to genotype-by-environment interactions, and the results here also reflect that it may not be feasible to set breeding for DM yield broad adaptation as a principle breeding goal. The comparatively low predictive abilities of the HA prediction models [20] applied (see below) are also likely to have influenced the overall selection responses observed for HA. This was likely affected by both the small size and the composite nature of the training set [23]. Future optimization of genomic selection should focus on increasing the training set size and a careful consideration of its composition [22] with regards to relatedness to the intended selection population(s). Expansion of the size of training sets for DM yield should be possible, with the advent of non-destructive sensorbased tools for measuring DM yield in plots [47,48], set to overcome phenotyping bottlenecks.
Our results, for HA as well as DTH, also indicated an overall difference in the selection response between Pop I and Pop III, with Pop I showing larger differentials between the Low and High GEBV SGs in both G1 and G3 ( Figure 3, Table 2). This suggests an influence of the genetic background on the GS response, although the underlying basis for this is not readily explained. Both Pop I and Pop III were equally represented in the original training set. The predictive abilities of the five HA prediction models, estimated by cross-validation in the training set [20], were not markedly different for the two populations, ranging from r = 0.08 to 0.30 in Pop I (mean r = 0.18) and r = 0.05 to 0.27 (mean r = 0.16) in Pop III. Although we assume that the predictive ability in the current study is principally underpinned by the exploitation of genomic relatedness amongst individuals [20], the extent of LD in the respective populations might also contribute to the selection response differences observed between Pop I and Pop III. Global LD was estimated to be considerably lower in Pop III (366 bp at r 2 = 0.25) compared to Pop I (1506 bp) in the training set generation [20], and this trend may be expected to have carried through to the generations used for GS, as indicated by the mean relatedness estimates (Table 1).

Conclusions
This study shows that GS can be effective in perennial ryegrass to improve traits with contrasting genetic architectures, including dry matter yield. The variation in the observed response to genomic selection highlights the importance of target environment selection for training prediction models, the influence of relatedness between the training set and selection populations, and a contribution of genetic background. Indicative GS responses observed in the G3 selection populations, which were generationally more distant to the training set, suggests that genomic selection may be effective over multiple cycles of recurrent selection, even when a small, composite population of plants is used to train the genomic prediction models. However, this, and the underlying causes for asymmetry and population-specificity in genomic selection responses, require further investigation.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1: Figure S1: Schematic showing the relationship between the perennial ryegrass training and selection populations used in this study. G1 = offpsring generation of the training set; G3 = advanced generation resulting from two cycles of conventional selection (for non-DM yield traits) and recombination. Training set populations Pop II, IV, and V were not subject to genomic selection in this study; genetic relationships amongst training set populations are described in [20]; Figure S2: Days-to-heading for Low GEBV, High GEBV, and Non-selected selection group (SG) progeny samples (n = 25 per SG) generated by genomic selection from four populations (Pop I G1, Pop I G3, Pop III G1, and Pop  III G3). Data are BLUEs across two seasons, for the performance at the two locations (a) Waikato and (b) Manawatu. Within each population, SG bars with different letters are significantly different (p < 0.05) based on least significant difference (L.S.D 0.05); Figure S3: Herbage accumulation (from n = 8 seasonal harvests across two environments) measured as the dry matter yield per plot, for Low GEBV, High GEBV, and Non-selected selection group (SG) progenies, generated by genomic selection from four populations (Pop I G1, Pop I G3, Pop III G1, and Pop III G3), and combined across all populations (All Pop). Data are BLUEs for the average performance across two locations (Waiakto and Manawatu). Within each population, SG bars with different letters are significantly different (p < 0.05), as supported by least significant difference (L.S.D 0.05) error bars; Table  S1: Schedule of measurements undertaken for herbage accumulation (HA, recorded as the dry matter yield per 2 m row) and days-to-heading (DTH, recorded as the number of days after 22 October) in trials located at Waikato and Manawatu; Table S2: Genomic relationship matrix.csv; Table S3: Genotypic (σ 2 g), genotype-byharvest (σ 2 gh), genotype-by-location (σ 2 gl), and experimental error (σ 2 ε) variance components and their associated standard errors (± SE) when analysed based on the selection group (Low GEBV, Non-selected, and High GEBV × four selection populations) and by selection category (selection groups combined across all populations) for mean herbage accumulation within-and across-location. Italicized numbers indicate a non-significant (p > 0.05) variance component.