BACKGROUND
perennial ryegrass  is one of the most cultivated forage species in temperate grasslands, mainly farmed for its re-growth capacity after defoliation, and for the high value as feed for ruminants, due to palatability, digestibility, and nutritive contents .
perennial ryegrass is an obligate allogamous species with genetic gametophytic self-incompatibility  <cit> , and is bred in genetically heterogeneous families.
heading date  is an important trait for forage species, often used as a model trait  <cit>  due to its high heritability and the ease of assess.
it follows the shift from vegetative to reproductive growth, and it is significantly correlated with several other traits involved in plant growth and development, such as plant height, spike length, tiller number and size, leaf length  <cit> , as well as with a number of yield and quality traits.
early heading genotypes show a higher growth rate in spring and higher forage yield in the first cut, .
correspondingly, the opposite was found for later cuts in summer, where dry matter forage yield was higher for late genotypes.
results in the literature are inconsistent regarding performances in fall and winter.
humphreys  <cit>  found higher autumn and winter growth rate in late genotypes, while in sampoux et al.
<cit> , the correlation between hd and forage yield in autumn was not significantly different from zero.
differences between early and late genotypes were also observed in the intensity of aftermath heading, which was higher in the early material  <cit> , as well as in the content of fiber and soluble sugars.
humphreys  <cit>  and sampoux et al.
<cit>  measured less water soluble carbohydrates and more lignin and neutral detergent fiber in early- than in late hd varieties.
late heading was also associated with higher digestibility and therefore to a higher lactation energy content for milk production  <cit> .
although heading also marks the production of seeds, the correlation with seed yield is unclear.
later genotypes were generally found to give lower seed production  <cit> , but this correlation was not always significant  <cit> .
due to its significant effects on other traits, breeding has always aimed to exploit the natural variation in hd, in order to create mixtures of varieties that could give high performances throughout the whole year.
since international listing of new varieties requires fulfilment of the three criteria; distinctiveness, uniformity, and stability  there is also a strong breeding focus on hd in order to create uniform varieties.
while the uniformity of inbred varieties is rather easy to control it can be more challenging in outbreeding grass varieties that are breed as families.
in order to ensure more stability in forage quality over the season, cultivars have been divided into different earliness groups.
the number- and extend of each hd group differs between countries, with some countries defining up to nine hd groups.
however, hd appears to behave as a continuous character, and the distinction between early and late material is not always clear, with new candidates that may be classified in different neighboring hd groups, depending on the definitions used in the different countries.
the trait always showed medium to high heritability .
kearsey et al.
<cit>  showed the presence of both additive and dominance effects, with the first being the larger and dominance being for early heading, but did not find any evidence for epistasis.
genotype by environment  interactions were found to be small by ravel and charmet  <cit> , in a multi-site analysis in france.
however, a different result was obtained by kearsey et al.
<cit> , who showed interactions between the environment and both additive and dominance effects, in an experiment across italy and the uk.
in the latter decade, the genetic control of hd was better understood thanks to the use of molecular markers and comparative genome analyzes.
in model species, such as arabidopsis thaliana l., as well as in cereals like wheat and rice heading or the control of flowering has been the subject for numerous studies and publications .
especially the use of induced arabidopsis mutants and the combinations of such lead to the detailed modelling of the genetic control of flowering in plants.
the investigations demonstrated the involvement of genes belonging to three major pathways:  vernalization response genes , which regulates heading after low temperature periods;  photoperiod response genes , which is active/inactive with a certain day length;  ‘earliness per se’ factors, which seems to be independent of light and cold requirements  <cit> .
in perennial ryegrass a number of flowering genes were previously identified by sequence homology with flowering genes found in arabidopsis, rice, and maize .
others were identified through classical quantitative trait loci  mapping, performed using different plant material and different genetic maps.
genetic maps were organized in seven linkage groups , numbered according to the conserved synteny with the triticeae’s maps  <cit> .
qtls were identified on all seven lgs .
comparison between studies is complicated due to lack of common markers and it is always difficult to determine if two significant markers found on the same lg, actually correspond to the same qtl.
furthermore, among different studies there is often poor agreement regarding the number and the distribution of the qtls, likely due to environmental factors, use of different mapping populations  <cit> , and low statistical power in several studies.
a great effort was put in understanding the genetic mechanisms behind the qtls in lg <dig> and lg <dig>  which were significant in almost all studies.
the first was found to be in a syntenic association with the wheat vrn <dig> gene  <cit> , and its function seems to be conserved between diploid wheat and perennial ryegrass  <cit> .
a relation was also hypothesized with a putative casein kinase gene, previously mapped in rice and involved in photoperiod sensitivity  <cit> .
the qtl on lg <dig> was suggested to be associated with the gene lpco, homologous to the constans of arabidopsis and the hd <dig> of rice, involved in the photoperiodic regulation of flowering time  <cit> .
synteny was also detected with the hd <dig> region of rice  <cit> , which codes for a flowering-locus-t  orthologue of arabidopsis.
ft gene is involved in induction to reproductive growth at the meristem  <cit>  and has been shown to actively regulate the flowering response in l. perenne  <cit> .
other hypothesis have been proposed to relate the other significant markers to qtl previously found in related species, such as lolium multiforum lam.
and festuca pratensis l.  <cit> .
while these studies identified some of the key genes in floral control in ryegrass, the biology of the trait is still far from being understood.
furthermore, the use of qtl analyses was shown to be not effective in capturing small effect genes  <cit>  and to overestimate the variance explained by qtls, due to the so called beavis effect  <cit> .
however, such limitations may be overcome by the use of genomic selection .
in contrast to traditional marker assisted selection , gs does not focus on finding specific qtls, but selects families/individuals based on genomic estimated breeding values , which are calculated using all markers simultaneously.
linkage disequilibrium  between causative loci and markers is ensured by high marker coverage.
such ld can come from three sources:  close physical linkage between marker and qtl;  family structure in the population, creating both short range  and long range  ld;  population structure due to mixing breeding material of different origin.
therefore, the ld can be also tracked across families, enabling to estimate marker effects at a population level  <cit> .
gs is practically implemented trough different steps:  model development on a set of individuals/families that are both genotyped and phenotyped ;  estimation of gebvs for a set of individuals/families that are only genotyped , based on their relationship with the training set;  selection of the best breeding material.
in this paper we will refer to the first two steps as genomic prediction .
gs is now widely used in animal breeding  <cit> , but it is still a new technology in crop breeding.
to date only a limited number of studies has been published on real data, mainly on species that are primarily grown as homogeneous varieties, such as maize, barley, and wheat .
the first results are promising and gs is expected to significantly increase genetic gains, especially due to the shortening of the breeding cycles  <cit> .
so far, aside from a few theoretical discussions, very little has been reported about gs potentials in allogamous species that are breed and farmed as heterogeneous populations.
specifically for perennial ryegrasses, hayes et al.
<cit>  showed good perspectives for introducing gs in practical breeding programs.
however, a full implementation would require radical changes in the present breeding systems, and may face problems due to low ld and high effective population size, due to the outcrossing nature of the species  <cit> .
this paper represents our first attempt to introduce gp in a breeding program of forage perennial ryegrass, using hd as model trait.
1757 f <dig> families , phenotyped for hd and genotyped with high marker coverage, were used to dissect the genetic and genomic structure of the trait.
first, a genome wide association analysis , to check for the presence of major qtl was conducted.
second, significant markers were used to calculate the gebvs in a set of synthetic  families, a part of which was related with the training set.
third, gp models were tested within the f <dig> set, using different cross-validation  schemes and different population sizes, and then used to predict breeding values of syn families.
predictive ability of gp was compared with predictions based on gwas results.
RESULTS
population structure, ld, and genetic parameters
results from the principal component analysis  showed the presence of some degree of population structure .
the ‘elbow’ point of the pca scree plot was determined at the fourth pc .
the first four pcs explained  <dig>   <dig>   <dig>  and 6% of the variance among snps respectively.
the optimal numbers of cluster, determined by k-means clustering, turned out to be two.
the separation in the two clusters could be explained by the first pc and it was strongly related to the origin of the parent populations .
in fig.
1a, all the families represented by blue points were identified as pair-crosses having a varieties originated in uk as one pp.
for this reason, in the following part of the paper, we will refer to this group as .
the population structure was also shown to be related to hd, which was mostly explained by the third pc .fig.
<dig> first pcs and: a origin of the pps; b breeding value for hd
results from the ld analysis are displayed in fig.
<dig>  the ld was shown to have a rather short extent, decaying below  <dig>  after a few hundred base pairs .
without using any correction for relatedness and population structure , the markers with ld >  <dig>  and ld >  <dig>  were the  <dig>  % and the  <dig>  % of the total number of snps respectively.
the average distance between markers having ld >  <dig>  was about 8900 bp, and for markers having ld >  <dig>  the distance was close to 3600 bp.
the correction further reduced the proportion of snps in ld , which dropped down to  <dig>  % for ld >  <dig>  and to  <dig>  % for ld >  <dig> , corresponding to a reduction of  <dig> and 40 % respectively.
the average distance between markers having ld >  <dig>  was reduced to  <dig>  while the one for markers having ld >  <dig> , drop down to 1200 bp, corresponding to a reduction of  <dig> and 66 % respectively, showing that the correction for structure and relatedness reduced the short range ld in the population.
the proportions of snps separated by more than  <dig> and 6300 bp were about  <dig> and 71 % respectively.fig.
<dig> ld decay: a without corrections; b corrected for relatedness and population structure
the total amount of phenotypic variance, together with different definition of heritability, is shown in fig.
<dig>  the additive genetic variance accounted for half of the total phenotypic variance, and it was equally divided between the ‘within pps’ and the ‘among pps’ components.
the interaction between the additive effect and the environment was relatively small  and occurred only within pps.fig.
<dig> variance components, phenotypic variance, and heritabilities .
legend: i = σ
2
i; p = σ
2
p; ily = σ
2
ily; pply = σ
2
pply; e = σ
2
e; h
2
n = narrow sense heritability across environments; h
2
ne = narrow sense environment-specific heritability; h
2
be broad sense environment-specific heritability
genome wide association
using the bayesian information criterion, the optimal number of pcs for population structure correction was determined to be four, confirming the visual identification of the ‘elbow’ point.
the effect of the correction on the significance levels expressed as –log <dig> is clear from the qq-plots reported in fig.
<dig>  after selection for high ld within scaffold, the number of significant snps  was  <dig> using the t-test with bonferroni correction, and  <dig> using false discovery rate  .
snps are anchored to genomic scaffolds, which are not orientated or ordered with respect to a genetic map.
however, the draft assembly has been annotated with the aid of extensive transcriptome data and a number of genes have been predicted in the scaffolds harboring the significant snps .
a total of ten markers were found to be within the gene space,  <dig> of which were mapped in exon regions .
the allele substitution effects ranged from  <dig>  to  <dig>  days.
the percentage of additive variance across locations/years explained by each marker was between  <dig>  and  <dig>  % within the f <dig> families, and between  <dig>  and  <dig>  % in the syn families.
the sum of the variances explained by all significant markers corresponded to about  <dig>  % in the f2s and  <dig>  % in the syns.
the correlation between the marker effect in the two sets was positive .fig.
<dig> qq-plots, without  and with  correction for g-matrix and pcs
the snp 5059| <dig> is situated in a scaffold where the hd <dig> homolog of the lpco gene was also mapped.
the marker 2801| <dig> locates in a gene encoding for phytocrome c .
the scaffold  <dig> harbored two significant snps: 1379| <dig>  situated in a gene encoding for a pectate lyase  <dig>  and 1379| <dig>  located outside the gene space, but only 3500 bp away from the pectate lyase  <dig> coding region.
prediction of syns based on snp markers is shown in fig.
<dig>  using only the most significant snp, it was possible to get an accuracy of predictive ability of  <dig> .
adding more markers, initially improved the predictions.
accuracy was  <dig>  using all the markers that were declared significant with the bonferroni corrected t-test.
the highest correlation between real and estimated breeding values was reached by using all markers that passed the significance threshold in the fdr test.
the estimate was equal to  <dig> , meaning that all significant snps were able to predict the  <dig>  % of the genetic variance in the syn families.
keeping adding markers after that threshold did not lead to any improvement in the predictions, giving accuracies that kept fluctuating between  <dig>  and  <dig> .
concerning the bias of the estimates, using only a few markers clearly led to a downward bias in the variance of the predicted breeding values.
the underestimation was on the order of  <dig> when the first  <dig> snps were used.
adding more markers led to a significant reduction in bias.
at the point of maximum predictive ability, it was  <dig> , and it kept decreasing even if insignificant markers were added.fig.
<dig> accuracies  and bias  for prediction of syns, based on marker effects.
legend: the vertical lines indicate the two significant thresholds: bonferroni corrected t-test , and fdr
genomic prediction
for each cross-validation  scheme, predictive abilities, accuracies, and bias are displayed in table  <dig>  within f2s, predictive abilities were extremely high, reaching the value of  <dig>  in the absence of related families in the two sets , and  <dig>  when training and validation set contained related families, in the so called k-fold scheme .
the hotelling-williams test showed the two values to be significantly different .
for the k-fold scheme, the accuracy computed using the approximation in  gave an estimate that was higher than one .
however, the standard error  of this estimate was  <dig> , indicating that the actual accuracy could range between  <dig>  and  <dig> .
analyses on reduced training sets  showed that accuracies above  <dig>  were reachable, for both schemes, with training populations as larger than  <dig> families.
bias was shown to be always very low in the k-fold scheme.
the situation was different in the pp-fold cv, where the gebvs variance was generally underestimated, and where an increase in the population size resulted in a bias reduction: the regression coefficient  was  <dig>  using  <dig> families, and  <dig>  using the full dataset.
bias for population sizes below  <dig> families is not shown, because affected by very large se and not indicative of any trend.
accuracies within the set  were lower than the ones found on an equal number of randomly chosen f2s, especially in the pp-fold scheme.
bias was also slightly higher.
the cv within the other cluster gave more or less the same results as the cv within all f2s.table  <dig> population size and results  for all cv schemes
†
a
b
a
b
a
b
n
n
a
f2s - > syns ‡
b
†different letters indicate a significant difference between the two cv schemes  based on hotelling-williams test.
n indicates that the comparison does not apply, as models were based on different sets of data
‡using all snps that were declared significant after fdr test
fig.
<dig> gebv vs. corrected mean phenotypes: a within f2s ; b predicting syns from f2s.
legend: blue line = plot diagonal; red line = linear regression
fig.
<dig> accuracies  and bias  with different population sizes: k-fold  and pp-fold scheme
predictions across sets also worked well.
accuracy of predicting uk set from the other f <dig> families was slightly lower .
predictions for gebvs were better when the set  was used as training set.
in this case, accuracies were comparable to the ones obtained within all f2s, with a pp-fold scheme, and using a similar population size.
the bias level indicates that the gebvs variance was underestimated when the set  was used as training population, and slightly overestimated when the prediction was performed in the opposite direction.
predictive ability for syns  was similar to the ones within all f2s, and significantly different from the one obtained from gwas results .
the accuracy was  <dig> , 14 % higher than the in prediction based on the significant markers.
in this case, the linear regression of mean corrected phenotypes on gebvs indicated no bias in the gebvs variance.
discussion
population structure, ld, and genetic parameters
the population structure was mainly defined by the origin of the pps, which was correlated with the first principal component.
the majority of the f <dig> families were grouped in one big cluster.
this may lead to the hypothesis of a common european genetic pool.
this pool is likely to originate from the continuous and  free exchange of breeding material among the different breeders.
the parents of the set  may be an exception to that pool, and their genetic origin need to be further investigated.
the relation between pc <dig> and hd indicates the need to correct for population structure while performing gwas, in order to avoid false positives.
further variance analyses were performed by adding fixed regressions for the first  <dig>   <dig>   <dig>  and  <dig> pcs to the equation shown in formula .
result indicates a correlation of hd with the pc <dig>  but not with the other three main pcs.
when accounting for the first two pcs, the additive genomic variance across location was equal to the 98 % of the ones of the model without any pc.
when accounting for the first three pcs, the additive genetic variance left was the 89 %.
adding a regression for pc <dig> had a negligible effect.
the ld within scaffolds showed to decay rapidly, confirming the concerns expressed by hayes et al.
<cit> , who reported useful ld  to extend at best 1 kb.
however, in the present breeding material, the presence of relatedness and population structure generally increased the ld, bringing to an increase by the order of three in the average distance between those markers.
this fact suggests that population structure, which is known to be mostly responsible for the long range ld, also plays an important role in increasing the level of ld within scaffolds.
the correction also led to a decrease in the average distance between markers in ld, which was more pronounced for higher ld levels.
estimation of variance components confirm results obtained by fè et al.
<cit>  on a subset of the same data.
in this paper was also possible to calculate the heritability across environments, and to estimate the extent of g × e for additive and non-additive effects.
compared with other traits previously analyzed  <cit> , the proportion of genetic variance between pps was much higher.
the small level of g × e seems to confirm the results obtained by ravel and charmet  <cit> .
however, plants were cultivated only in denmark and england.
to have a better understanding of g × e effects, it would be a good idea to perform experiments covering more diverse climatic conditions.
the significant amount of σ2pply may indicate the presence of dominance acting between families and within single environments.
in the literature, the presence of non-additive effect is reported also across location  <cit> .
however results are not directly comparable, as this paper ignores additive effects that may be present within f <dig> families.
significant markers and genetic architecture
the gwas analysis revealed a rather complex genetic architecture of hd in ryegrass.
several markers with significant effect were identified.
in the most significant snps, a shift from one homozygous form to the other can cause changes of up to  <dig>  days in the date of heading.
that is a remarkable difference, if compared with the level of variation in the phenotypes: average phenotypes corrected for fixed effect had a sd of  <dig> .
however, due to low minor allele frequency , these markers were only able to explain a small proportion of the additive variance, which may indicate the presence of a large number of genes also affecting the trait, but with effects lower than the detection limit.
there were a high number of significant snps found outside the gene space.
this is not too surprising considering a recent study in maize found that the majority of trait associated variants were located outside annotated genes, but within  <dig> kb of transcriptional start and stop sites  <cit> .
some of the significant snps were clearly linked to genes that may have a direct or indirect influence on hd.
it is well established that constans  plays a crucial role in promoting flowering in response to long days  <cit> .
we identified a significant snp within less than  <dig> kb of a co homolog.
it has already been established that a ryegrass homolog to co exhibits expression patterns consistent with its function in arabidopsis, and can complement co mutants  <cit> .
furthermore, the co homolog co-located on linkage group seven with a large effect qtl for hd  <cit> .
allelic variation in an intergenic region upstream of co was found to be significantly associated with hd in a collection of  <dig> perennial ryegrass genotypes   <cit> .
the ~ <dig> kb region sequenced as part of that study shares near perfect identity with scaffold  <dig> , which has the significant snp identified in our study, and therefore represent the same genomic regions.
overall, our results provide further evidence that allelic variation at co contributes to variation in hd in perennial ryegrass, specifically within a large collection of breeding families.
we also identified a significant snp within the coding region of a homolog to phyc.
phytochromes are red/far-red photoreceptors that play a role in how a plant responds to light, and adapts its growth and development.
it was recently demonstrated in wheat that phyc plays a major role in accelerating flowering under long-days  <cit> , in contrast to the model plants such as arabidopsis and rice where phyc represses flowering under non-inductive conditions.
the fact that loss-of-function mutations in phyc resulted on average in a 108 day delay in flowering of wheat under long days emphasizes the potential for allelic variation at this gene to greatly alter flowering times.
a similar role for phyc has also been recently reported in brachypodium distachyon  <cit> .
perennial ryegrass is a close relative of wheat and brachypodium, and a similar role for phyc in floral induction of perennial ryegrass is possible.
a homolog to phyc has been mapped to linkage group four of perennial ryegrass  <cit> , although it mapped some distance from the hd qtl identified in that experimental population.
no correlation was found between significant markers and other genes that are known to be important in flowering time regulation, such as ft. that may be due to different reasons:  absence of causative polymorphisms in the breeding material;  no or low ld between markers and the causative polymorphisms ;  low maf at the causative polymorphisms ;  polymorphisms not detected because they are correlated with the family structure and shrunken by the correction with g-matrix and pcs .
prediction of breeding values
despite explaining only a small part of the genetic variance in the syns, the significant markers were able to predict the breeding values with high accuracy, even when only a few markers were used.
this is due to few genes with relatively large effects identified in the f <dig> population.
however, the presence of a certain level of population structure  will also contribute to the predictive ability in the syns.
in the gwas, we accounted for the presence of population structure by correcting the marker effect .
however, that correction does not apply to the estimation of prediction accuracy.
when we correlate the phenotypes with one marker, we are actually estimating the correlation of the phenotype with that particular marker, plus all the population structure that is correlated to the snp.
the trend in accuracy for an increasing number of markers met our expectations: any significant snp is supposed to add information that will increase the correlation with the true breeding value.
non-significant snps will mainly add random noise to the correlation, but were able to add genetic information that increased the variance of estimated breeding values.
the fact that the accuracy reaches the highest value in correspondence of the nineteenth snp is also a strong argument for using fdr, instead of bonferroni corrected t-test, as significance test.
the decreases in prediction accuracy that happened after adding the fourth and the sixth markers may be related to different levels of expression or to different interactions in the two populations.
results from syns prediction  show a clear advantage for using gp, compared with gwas, both in terms of accuracies, and in terms of bias, as well as its good potential in predicting across different generations.
the relatively high se for accuracies may be related to the experimental design that, due to incomplete randomization between trials and pps, which could lead to less accurate estimate of pp variance components.
within all f <dig> families, predictions were extremely good, allowing the explanation of nearly the whole genetic variance.
this result is higher than what usually found for the same trait in other species , even though accuracies above  <dig>  have already been reported in other outcrossing species such as maize.
that may be related with the high level of structure in the population, and to the fact that heading date primarily is affected by additive genetic effects, so the additive values of the pps are very well estimated.
a very high accuracy may also seem in contrast with what reported in the literature for traits affected by major effect snps  <cit> .
theoretically, for traits that include some genes with large effect, it would be recommended to use other prediction methods such as bayesian models, which allow marker effects to belong to distributions with different variance.
however, genomic best linear unbiased prediction , when compared with bayesian methods, was shown to be better in accounting for population structure, but less capable to explain the short range ld between markers  <cit> .
this makes it particularly effective for gp in breeding programs of species like perennial ryegrass, characterized quick decay of short range ld, and usually bred on a sib-mating scheme.
the lower accuracies found within the set  may be related to a low level of population structure within the cluster, as appear also from fig.
<dig>
accuracies reached by predicting from related families  were significantly higher than the ones obtained in the absence of related families , for any population size.
regarding the relationship between predictions and the size of training population, increasing the training size led to pronounced gains in predictive abilities for population sizes lower than  <dig> families, and to smaller gains in the case of larger populations.
despite that, due to the higher predictive abilities, accurate predictions could be obtained even using a relatively small training set.
problems of underestimations of the gebvs variance may arise in absence of closely related families.
regarding predictions across sets, the relative difficulty in predicting the set  may be due to families in this set having fewer relatives in the chosen training population.
the level of bias may arise due to the differences in genetic variance or to genetic correlations that are less than unity between the two sets.
this problem could compromise a correct comparisons between gebvs from the training set  and the validation set .
this should be taken into consideration during the design of the training population, allowing the presence a wide variety of genotypes in as many environments as possible.
CONCLUSIONS
our research clearly showed considerable potential for implementation of gs in breeding of l. perenne.
results obtained by gp significantly outperformed the accuracy based on traditional mas, being able to predict a very large proportion of the genetic variance.
gblup was shown to be capable of reaching very high accuracies, even in a trait characterized by major effect genes, at least in a population with fast decaying ld and population structure arising from admixture and relatedness.
predictions were also very good across datasets, with accuracies of up to  <dig> .
bias in the gebv variance could be caused by lack of common parent populations between training and validation set.
the study has also revealed important details about the genetic architecture of hd in l. perenne.
the trait appears to be controlled by both major and minor effect genes, regulated both by sequence changes within coding regions, and by the action of intergenic regulatory elements.
snps were identified within or proximal to genes with well-established roles in floral induction in plants .
despite this, the technique used for gwas has limitations, mainly due to the marker density given the rapid decay of ld, and due to the strong structure in the population.

 perennial ryegrass is one of the most cultivated forage species in temperate grasslands, mainly farmed for its re-growth capacity after defoliation, and for the high value as feed for ruminants, due to palatability, digestibility, and nutritive contents. perennial ryegrass is an obligate allogamous species with genetic gametophytic self-incompatibility, and is bred in genetically heterogeneous families. heading date is an important trait for forage species, often used as a model trait due to its high heritability and the ease of assess. it follows the shift from vegetative to reproductive growth, and it is significantly correlated with several other traits involved in plant growth and development, such as plant height, spike length, tiller number and size, leaf length, as well as with a number of yield and quality traits. early heading genotypes show a higher growth rate in spring and higher forage yield in the first cut,. correspondingly, the opposite was found for later cuts in summer, where dry matter forage yield was higher for late genotypes. results in the literature are inconsistent regarding performances in fall and winter. humphreys found higher autumn and winter growth rate in late genotypes, while in sampoux et al., the correlation between hd and forage yield in autumn was not significantly different from zero. differences between early and late genotypes were also observed in the intensity of aftermath heading, which was higher in the early material, as well as in the content of fiber and soluble sugars. humphreys and sampoux et al. measured less water soluble carbohydrates and more lignin and neutral detergent fiber in early- than in late hd varieties. late heading was also associated with higher digestibility and therefore to a higher lactation energy content for milk production. although heading also marks the production of seeds, the correlation with seed yield is unclear. later genotypes were generally found to give lower seed production, but this correlation was not always significant. due to its significant effects on other traits, breeding has always aimed to exploit the natural variation in hd, in order to create mixtures of varieties that could give high performances throughout the whole year. since international listing of new varieties requires fulfilment of the three criteria; distinctiveness, uniformity, and stability there is also a strong breeding focus on hd in order to create uniform varieties. while the uniformity of inbred varieties is rather easy to control it can be more challenging in outbreeding grass varieties that are breed as families. in order to ensure more stability in forage quality over the season, cultivars have been divided into different earliness groups. the number- and extend of each hd group differs between countries, with some countries defining up to nine hd groups. however, hd appears to behave as a continuous character, and the distinction between early and late material is not always clear, with new candidates that may be classified in different neighboring hd groups, depending on the definitions used in the different countries. the trait always showed medium to high heritability. kearsey et al. showed the presence of both additive and dominance effects, with the first being the larger and dominance being for early heading, but did not find any evidence for epistasis. genotype by environment interactions were found to be small by ravel and charmet, in a multi-site analysis in france. however, a different result was obtained by kearsey et al., who showed interactions between the environment and both additive and dominance effects, in an experiment across italy and the uk. in the latter decade, the genetic control of hd was better understood thanks to the use of molecular markers and comparative genome analyzes. in model species, such as arabidopsis thaliana l., as well as in cereals like wheat and rice heading or the control of flowering has been the subject for numerous studies and publications. especially the use of induced arabidopsis mutants and the combinations of such lead to the detailed modelling of the genetic control of flowering in plants. the investigations demonstrated the involvement of genes belonging to three major pathways: vernalization response genes, which regulates heading after low temperature periods; photoperiod response genes, which is active/inactive with a certain day length; ‘earliness per se’ factors, which seems to be independent of light and cold requirements. in perennial ryegrass a number of flowering genes were previously identified by sequence homology with flowering genes found in arabidopsis, rice, and maize. others were identified through classical quantitative trait loci mapping, performed using different plant material and different genetic maps. genetic maps were organized in seven linkage groups, numbered according to the conserved synteny with the triticeae’s maps. qtls were identified on all seven lgs. comparison between studies is complicated due to lack of common markers and it is always difficult to determine if two significant markers found on the same lg, actually correspond to the same qtl. furthermore, among different studies there is often poor agreement regarding the number and the distribution of the qtls, likely due to environmental factors, use of different mapping populations, and low statistical power in several studies. a great effort was put in understanding the genetic mechanisms behind the qtls in lg and lg which were significant in almost all studies. the first was found to be in a syntenic association with the wheat vrn gene, and its function seems to be conserved between diploid wheat and perennial ryegrass. a relation was also hypothesized with a putative casein kinase gene, previously mapped in rice and involved in photoperiod sensitivity. the qtl on lg was suggested to be associated with the gene lpco, homologous to the constans of arabidopsis and the hd of rice, involved in the photoperiodic regulation of flowering time. synteny was also detected with the hd region of rice, which codes for a flowering-locus-t orthologue of arabidopsis. ft gene is involved in induction to reproductive growth at the meristem and has been shown to actively regulate the flowering response in l. perenne. other hypothesis have been proposed to relate the other significant markers to qtl previously found in related species, such as lolium multiforum lam. and festuca pratensis l.. while these studies identified some of the key genes in floral control in ryegrass, the biology of the trait is still far from being understood. furthermore, the use of qtl analyses was shown to be not effective in capturing small effect genes and to overestimate the variance explained by qtls, due to the so called beavis effect. however, such limitations may be overcome by the use of genomic selection. in contrast to traditional marker assisted selection, gs does not focus on finding specific qtls, but selects families/individuals based on genomic estimated breeding values, which are calculated using all markers simultaneously. linkage disequilibrium between causative loci and markers is ensured by high marker coverage. such ld can come from three sources: close physical linkage between marker and qtl; family structure in the population, creating both short range and long range ld; population structure due to mixing breeding material of different origin. therefore, the ld can be also tracked across families, enabling to estimate marker effects at a population level. gs is practically implemented trough different steps: model development on a set of individuals/families that are both genotyped and phenotyped; estimation of gebvs for a set of individuals/families that are only genotyped, based on their relationship with the training set; selection of the best breeding material. in this paper we will refer to the first two steps as genomic prediction. gs is now widely used in animal breeding, but it is still a new technology in crop breeding. to date only a limited number of studies has been published on real data, mainly on species that are primarily grown as homogeneous varieties, such as maize, barley, and wheat. the first results are promising and gs is expected to significantly increase genetic gains, especially due to the shortening of the breeding cycles. so far, aside from a few theoretical discussions, very little has been reported about gs potentials in allogamous species that are breed and farmed as heterogeneous populations. specifically for perennial ryegrasses, hayes et al. showed good perspectives for introducing gs in practical breeding programs. however, a full implementation would require radical changes in the present breeding systems, and may face problems due to low ld and high effective population size, due to the outcrossing nature of the species. this paper represents our first attempt to introduce gp in a breeding program of forage perennial ryegrass, using hd as model trait. 1757 f families, phenotyped for hd and genotyped with high marker coverage, were used to dissect the genetic and genomic structure of the trait. first, a genome wide association analysis, to check for the presence of major qtl was conducted. second, significant markers were used to calculate the gebvs in a set of synthetic families, a part of which was related with the training set. third, gp models were tested within the f set, using different cross-validation schemes and different population sizes, and then used to predict breeding values of syn families. predictive ability of gp was compared with predictions based on gwas results.  population structure, ld, and genetic parameters results from the principal component analysis showed the presence of some degree of population structure. the ‘elbow’ point of the pca scree plot was determined at the fourth pc. the first four pcs explained and 6% of the variance among snps respectively. the optimal numbers of cluster, determined by k-means clustering, turned out to be two. the separation in the two clusters could be explained by the first pc and it was strongly related to the origin of the parent populations. in fig. 1a, all the families represented by blue points were identified as pair-crosses having a varieties originated in uk as one pp. for this reason, in the following part of the paper, we will refer to this group as. the population structure was also shown to be related to hd, which was mostly explained by the third pc.fig. first pcs and: a origin of the pps; b breeding value for hd results from the ld analysis are displayed in fig.  the ld was shown to have a rather short extent, decaying below after a few hundred base pairs. without using any correction for relatedness and population structure, the markers with ld >  and ld >  were the  % and the  % of the total number of snps respectively. the average distance between markers having ld >  was about 8900 bp, and for markers having ld >  the distance was close to 3600 bp. the correction further reduced the proportion of snps in ld, which dropped down to  % for ld >  and to  % for ld > , corresponding to a reduction of and 40 % respectively. the average distance between markers having ld >  was reduced to while the one for markers having ld > , drop down to 1200 bp, corresponding to a reduction of and 66 % respectively, showing that the correction for structure and relatedness reduced the short range ld in the population. the proportions of snps separated by more than and 6300 bp were about and 71 % respectively.fig. ld decay: a without corrections; b corrected for relatedness and population structure the total amount of phenotypic variance, together with different definition of heritability, is shown in fig.  the additive genetic variance accounted for half of the total phenotypic variance, and it was equally divided between the ‘within pps’ and the ‘among pps’ components. the interaction between the additive effect and the environment was relatively small and occurred only within pps.fig. variance components, phenotypic variance, and heritabilities. legend: i = σ 2 i; p = σ 2 p; ily = σ 2 ily; pply = σ 2 pply; e = σ 2 e; h 2 n = narrow sense heritability across environments; h 2 ne = narrow sense environment-specific heritability; h 2 be broad sense environment-specific heritability genome wide association using the bayesian information criterion, the optimal number of pcs for population structure correction was determined to be four, confirming the visual identification of the ‘elbow’ point. the effect of the correction on the significance levels expressed as –log is clear from the qq-plots reported in fig.  after selection for high ld within scaffold, the number of significant snps was using the t-test with bonferroni correction, and using false discovery rate. snps are anchored to genomic scaffolds, which are not orientated or ordered with respect to a genetic map. however, the draft assembly has been annotated with the aid of extensive transcriptome data and a number of genes have been predicted in the scaffolds harboring the significant snps. a total of ten markers were found to be within the gene space, of which were mapped in exon regions. the allele substitution effects ranged from to  days. the percentage of additive variance across locations/years explained by each marker was between and  % within the f families, and between and  % in the syn families. the sum of the variances explained by all significant markers corresponded to about  % in the f2s and  % in the syns. the correlation between the marker effect in the two sets was positive.fig. qq-plots, without and with correction for g-matrix and pcs the snp 5059| is situated in a scaffold where the hd homolog of the lpco gene was also mapped. the marker 2801| locates in a gene encoding for phytocrome c. the scaffold harbored two significant snps: 1379| situated in a gene encoding for a pectate lyase and 1379| located outside the gene space, but only 3500 bp away from the pectate lyase coding region. prediction of syns based on snp markers is shown in fig.  using only the most significant snp, it was possible to get an accuracy of predictive ability of. adding more markers, initially improved the predictions. accuracy was using all the markers that were declared significant with the bonferroni corrected t-test. the highest correlation between real and estimated breeding values was reached by using all markers that passed the significance threshold in the fdr test. the estimate was equal to, meaning that all significant snps were able to predict the  % of the genetic variance in the syn families. keeping adding markers after that threshold did not lead to any improvement in the predictions, giving accuracies that kept fluctuating between and. concerning the bias of the estimates, using only a few markers clearly led to a downward bias in the variance of the predicted breeding values. the underestimation was on the order of when the first snps were used. adding more markers led to a significant reduction in bias. at the point of maximum predictive ability, it was, and it kept decreasing even if insignificant markers were added.fig. accuracies and bias for prediction of syns, based on marker effects. legend: the vertical lines indicate the two significant thresholds: bonferroni corrected t-test, and fdr genomic prediction for each cross-validation scheme, predictive abilities, accuracies, and bias are displayed in table  within f2s, predictive abilities were extremely high, reaching the value of in the absence of related families in the two sets, and when training and validation set contained related families, in the so called k-fold scheme. the hotelling-williams test showed the two values to be significantly different. for the k-fold scheme, the accuracy computed using the approximation in gave an estimate that was higher than one. however, the standard error of this estimate was, indicating that the actual accuracy could range between and. analyses on reduced training sets showed that accuracies above were reachable, for both schemes, with training populations as larger than families. bias was shown to be always very low in the k-fold scheme. the situation was different in the pp-fold cv, where the gebvs variance was generally underestimated, and where an increase in the population size resulted in a bias reduction: the regression coefficient was using families, and using the full dataset. bias for population sizes below families is not shown, because affected by very large se and not indicative of any trend. accuracies within the set were lower than the ones found on an equal number of randomly chosen f2s, especially in the pp-fold scheme. bias was also slightly higher. the cv within the other cluster gave more or less the same results as the cv within all f2s.table population size and results for all cv schemes † a b a b a b n n a f2s - > syns ‡ b †different letters indicate a significant difference between the two cv schemes based on hotelling-williams test. n indicates that the comparison does not apply, as models were based on different sets of data ‡using all snps that were declared significant after fdr test fig. gebv vs. corrected mean phenotypes: a within f2s; b predicting syns from f2s. legend: blue line = plot diagonal; red line = linear regression fig. accuracies and bias with different population sizes: k-fold and pp-fold scheme predictions across sets also worked well. accuracy of predicting uk set from the other f families was slightly lower. predictions for gebvs were better when the set was used as training set. in this case, accuracies were comparable to the ones obtained within all f2s, with a pp-fold scheme, and using a similar population size. the bias level indicates that the gebvs variance was underestimated when the set was used as training population, and slightly overestimated when the prediction was performed in the opposite direction. predictive ability for syns was similar to the ones within all f2s, and significantly different from the one obtained from gwas results. the accuracy was, 14 % higher than the in prediction based on the significant markers. in this case, the linear regression of mean corrected phenotypes on gebvs indicated no bias in the gebvs variance. discussion population structure, ld, and genetic parameters the population structure was mainly defined by the origin of the pps, which was correlated with the first principal component. the majority of the f families were grouped in one big cluster. this may lead to the hypothesis of a common european genetic pool. this pool is likely to originate from the continuous and free exchange of breeding material among the different breeders. the parents of the set may be an exception to that pool, and their genetic origin need to be further investigated. the relation between pc and hd indicates the need to correct for population structure while performing gwas, in order to avoid false positives. further variance analyses were performed by adding fixed regressions for the first and pcs to the equation shown in formula. result indicates a correlation of hd with the pc but not with the other three main pcs. when accounting for the first two pcs, the additive genomic variance across location was equal to the 98 % of the ones of the model without any pc. when accounting for the first three pcs, the additive genetic variance left was the 89 %. adding a regression for pc had a negligible effect. the ld within scaffolds showed to decay rapidly, confirming the concerns expressed by hayes et al., who reported useful ld to extend at best 1 kb. however, in the present breeding material, the presence of relatedness and population structure generally increased the ld, bringing to an increase by the order of three in the average distance between those markers. this fact suggests that population structure, which is known to be mostly responsible for the long range ld, also plays an important role in increasing the level of ld within scaffolds. the correction also led to a decrease in the average distance between markers in ld, which was more pronounced for higher ld levels. estimation of variance components confirm results obtained by fè et al. on a subset of the same data. in this paper was also possible to calculate the heritability across environments, and to estimate the extent of g × e for additive and non-additive effects. compared with other traits previously analyzed, the proportion of genetic variance between pps was much higher. the small level of g × e seems to confirm the results obtained by ravel and charmet. however, plants were cultivated only in denmark and england. to have a better understanding of g × e effects, it would be a good idea to perform experiments covering more diverse climatic conditions. the significant amount of σ2pply may indicate the presence of dominance acting between families and within single environments. in the literature, the presence of non-additive effect is reported also across location. however results are not directly comparable, as this paper ignores additive effects that may be present within f families. significant markers and genetic architecture the gwas analysis revealed a rather complex genetic architecture of hd in ryegrass. several markers with significant effect were identified. in the most significant snps, a shift from one homozygous form to the other can cause changes of up to  days in the date of heading. that is a remarkable difference, if compared with the level of variation in the phenotypes: average phenotypes corrected for fixed effect had a sd of. however, due to low minor allele frequency, these markers were only able to explain a small proportion of the additive variance, which may indicate the presence of a large number of genes also affecting the trait, but with effects lower than the detection limit. there were a high number of significant snps found outside the gene space. this is not too surprising considering a recent study in maize found that the majority of trait associated variants were located outside annotated genes, but within kb of transcriptional start and stop sites. some of the significant snps were clearly linked to genes that may have a direct or indirect influence on hd. it is well established that constans plays a crucial role in promoting flowering in response to long days. we identified a significant snp within less than kb of a co homolog. it has already been established that a ryegrass homolog to co exhibits expression patterns consistent with its function in arabidopsis, and can complement co mutants. furthermore, the co homolog co-located on linkage group seven with a large effect qtl for hd. allelic variation in an intergenic region upstream of co was found to be significantly associated with hd in a collection of perennial ryegrass genotypes. the ~ kb region sequenced as part of that study shares near perfect identity with scaffold, which has the significant snp identified in our study, and therefore represent the same genomic regions. overall, our results provide further evidence that allelic variation at co contributes to variation in hd in perennial ryegrass, specifically within a large collection of breeding families. we also identified a significant snp within the coding region of a homolog to phyc. phytochromes are red/far-red photoreceptors that play a role in how a plant responds to light, and adapts its growth and development. it was recently demonstrated in wheat that phyc plays a major role in accelerating flowering under long-days, in contrast to the model plants such as arabidopsis and rice where phyc represses flowering under non-inductive conditions. the fact that loss-of-function mutations in phyc resulted on average in a 108 day delay in flowering of wheat under long days emphasizes the potential for allelic variation at this gene to greatly alter flowering times. a similar role for phyc has also been recently reported in brachypodium distachyon. perennial ryegrass is a close relative of wheat and brachypodium, and a similar role for phyc in floral induction of perennial ryegrass is possible. a homolog to phyc has been mapped to linkage group four of perennial ryegrass, although it mapped some distance from the hd qtl identified in that experimental population. no correlation was found between significant markers and other genes that are known to be important in flowering time regulation, such as ft. that may be due to different reasons: absence of causative polymorphisms in the breeding material; no or low ld between markers and the causative polymorphisms; low maf at the causative polymorphisms; polymorphisms not detected because they are correlated with the family structure and shrunken by the correction with g-matrix and pcs. prediction of breeding values despite explaining only a small part of the genetic variance in the syns, the significant markers were able to predict the breeding values with high accuracy, even when only a few markers were used. this is due to few genes with relatively large effects identified in the f population. however, the presence of a certain level of population structure will also contribute to the predictive ability in the syns. in the gwas, we accounted for the presence of population structure by correcting the marker effect. however, that correction does not apply to the estimation of prediction accuracy. when we correlate the phenotypes with one marker, we are actually estimating the correlation of the phenotype with that particular marker, plus all the population structure that is correlated to the snp. the trend in accuracy for an increasing number of markers met our expectations: any significant snp is supposed to add information that will increase the correlation with the true breeding value. non-significant snps will mainly add random noise to the correlation, but were able to add genetic information that increased the variance of estimated breeding values. the fact that the accuracy reaches the highest value in correspondence of the nineteenth snp is also a strong argument for using fdr, instead of bonferroni corrected t-test, as significance test. the decreases in prediction accuracy that happened after adding the fourth and the sixth markers may be related to different levels of expression or to different interactions in the two populations. results from syns prediction show a clear advantage for using gp, compared with gwas, both in terms of accuracies, and in terms of bias, as well as its good potential in predicting across different generations. the relatively high se for accuracies may be related to the experimental design that, due to incomplete randomization between trials and pps, which could lead to less accurate estimate of pp variance components. within all f families, predictions were extremely good, allowing the explanation of nearly the whole genetic variance. this result is higher than what usually found for the same trait in other species, even though accuracies above have already been reported in other outcrossing species such as maize. that may be related with the high level of structure in the population, and to the fact that heading date primarily is affected by additive genetic effects, so the additive values of the pps are very well estimated. a very high accuracy may also seem in contrast with what reported in the literature for traits affected by major effect snps. theoretically, for traits that include some genes with large effect, it would be recommended to use other prediction methods such as bayesian models, which allow marker effects to belong to distributions with different variance. however, genomic best linear unbiased prediction, when compared with bayesian methods, was shown to be better in accounting for population structure, but less capable to explain the short range ld between markers. this makes it particularly effective for gp in breeding programs of species like perennial ryegrass, characterized quick decay of short range ld, and usually bred on a sib-mating scheme. the lower accuracies found within the set may be related to a low level of population structure within the cluster, as appear also from fig.  accuracies reached by predicting from related families were significantly higher than the ones obtained in the absence of related families, for any population size. regarding the relationship between predictions and the size of training population, increasing the training size led to pronounced gains in predictive abilities for population sizes lower than families, and to smaller gains in the case of larger populations. despite that, due to the higher predictive abilities, accurate predictions could be obtained even using a relatively small training set. problems of underestimations of the gebvs variance may arise in absence of closely related families. regarding predictions across sets, the relative difficulty in predicting the set may be due to families in this set having fewer relatives in the chosen training population. the level of bias may arise due to the differences in genetic variance or to genetic correlations that are less than unity between the two sets. this problem could compromise a correct comparisons between gebvs from the training set and the validation set. this should be taken into consideration during the design of the training population, allowing the presence a wide variety of genotypes in as many environments as possible.  our research clearly showed considerable potential for implementation of gs in breeding of l. perenne. results obtained by gp significantly outperformed the accuracy based on traditional mas, being able to predict a very large proportion of the genetic variance. gblup was shown to be capable of reaching very high accuracies, even in a trait characterized by major effect genes, at least in a population with fast decaying ld and population structure arising from admixture and relatedness. predictions were also very good across datasets, with accuracies of up to. bias in the gebv variance could be caused by lack of common parent populations between training and validation set. the study has also revealed important details about the genetic architecture of hd in l. perenne. the trait appears to be controlled by both major and minor effect genes, regulated both by sequence changes within coding regions, and by the action of intergenic regulatory elements. snps were identified within or proximal to genes with well-established roles in floral induction in plants. despite this, the technique used for gwas has limitations, mainly due to the marker density given the rapid decay of ld, and due to the strong structure in the population. 
then, genomic prediction models were tested and prediction accuracy was compared to the one obtained in traditional marker assisted selection methods. we have demonstrated that there are good perspectives for gs implementation in perennial ryegrass breeding, and that problems resulting from low linkage disequilibrium can be reduced by the presence of structure and related families in the breeding population. genome wide association analysis was performed to uncover the genetic architecture of the trait. while comprehensive genome wide association analysis is difficult in species with extremely low ld, we did identify variants proximal to genes with a known role in flowering time. here, we test the opportunity to implement gs in breeding of perennial ryegrass, using real data from a forage breeding program. genomic selectionperennial ryegrassheadingfloweringgwasloliumconstansphytochromeaccuracyissue-copyright-statement© the author 2015 genomic selection has become a commonly used technology in animal breeding. so far, its implementation in plant breeding has been mainly investigated in species farmed as homogeneous varieties. gp models gave very high accuracies, which were significantly better than those obtained through traditional mas. accuracies were higher when predictions were made from related families and from larger training populations, whereas predicting from unrelated families caused the variance of the estimated breeding values to be biased downwards. heading date was chosen as a model trait, due to its high heritability and ease of assessment.
 genomic selection has become a commonly used technology in animal breeding. in crops, it is expected to significantly improve the genetic gains per unit of time. so far, its implementation in plant breeding has been mainly investigated in species farmed as homogeneous varieties. concerning crops farmed in family pools, only a few theoretical studies are currently available. here, we test the opportunity to implement gs in breeding of perennial ryegrass, using real data from a forage breeding program. heading date was chosen as a model trait, due to its high heritability and ease of assessment. genome wide association analysis was performed to uncover the genetic architecture of the trait. then, genomic prediction models were tested and prediction accuracy was compared to the one obtained in traditional marker assisted selection methods.  several markers were significantly associated with heading date, some locating within or proximal to genes with a well-established role in floral regulation. gp models gave very high accuracies, which were significantly better than those obtained through traditional mas. accuracies were higher when predictions were made from related families and from larger training populations, whereas predicting from unrelated families caused the variance of the estimated breeding values to be biased downwards.  we have demonstrated that there are good perspectives for gs implementation in perennial ryegrass breeding, and that problems resulting from low linkage disequilibrium can be reduced by the presence of structure and related families in the breeding population. while comprehensive genome wide association analysis is difficult in species with extremely low ld, we did identify variants proximal to genes with a known role in flowering time. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords genomic selectionperennial ryegrassheadingfloweringgwasloliumconstansphytochromeaccuracyissue-copyright-statement© the author 2015 