Phenotypic diversity

Since grasslands need to be adapted to the changing climate while maintaining a sufficient agronomical value, we aimed to evaluate the natural genetic diversity of perennial ryegrass for traits expected to contribute to environmental adaptation and/or agronomical value. Based on literature about environmental adaptation of grass species, some traits could be presumed as adaptive. Those include autumn growth rate, reproductive tiller density, heading earliness (Barre et al. 2018), isotopic composition of ¹³C (Condon et al. 2002), vegetative plant height before stem elongation (Oyarzabal et al. 2007), aerial dry matter lignin content (Casler et al. 2002), water soluble carbohydrates content of lamina or aerial dry matter (Humphreys 1989; Volaire 1995; Sanada et al. 2007). In our study, different patterns of spatial autocorrelation were observed according to traits (Table S7). Those spatial autocorrelations represent genetic patterns since they were calculated using the coordinates of the sites of origin of the populations and phenotypic data recorded at common trial sites. This suggests that demographic and selection processes have not similarly impacted the genetic diversity involved in the determinism of these various traits. As it has been proposed that most of the natural phenotypic diversity is due to more or less ancient Darwinian selection (Rieseberg et al. 2002), most traits for which wide natural genetic diversity has been observed may be adaptive in some way.

To apprehend the services that raw natural diversity can render to farmers, traits of interest for forage production were evaluated. Those include forage quality and phenological traits, disease resistance, vigor after cutting and growth at various seasons assessed through the measurement of canopy heights. We found that canopy height was fairly well correlated (File S1) to dry matter yield (r > 0.8) recorded at the same date during the spring vegetative growth, which confirms that canopy height could be used as an indicator of biomass production (Powell 1974).

The relatively high H² values observed for most traits (Table S7) indicated that great inter-population variance was revealed by the field trials implemented for this study. For example, heading date heritability was greater than 0.97 in all environments, which is higher than values reported in previous studies (Fè et al. 2015b; Barre et al. 2018). The dry matter yields measured on a sample of plots at each site (File S1) were within a common range for perennial ryegrass (McDonagh et al. 2016).

Population × environment interactions were highly significant for all recorded traits, even those linked to phenology such as heading date (Table S7). This could be explained by the marked differences of climate between the three experimental sites. For example, the spring growth (weekly records of canopy height) in the exceptionally cold conditions at PO in 2016 (Table S5) was poorly correlated to that observed in other locations or in other spring season at PO (Table S8). Adaptation to frost stress was likely revealed for some populations that grew best in spring 2016 at PO since they also suffered the less damage in the preceding winter (2015-16). A previous study which evaluated perennial ryegrass natural populations in experimental sites located in France also revealed strong population × environment interactions even for phenological traits such as heading date (Charmet et al. 1990).

Traits important for forage production often displayed different phenotypic distributions between the set of cultivars and the set of natural populations (Table S7). The values of cultivars were often tightly distributed closer toward an interesting extreme for forage production. However, for all traits, the range of values shown by the cultivars was within (or close to) the range of values shown by the natural populations. This was even the case for productivity traits such as spring canopy heights or autumn growth rates and for forage quality traits such as lignin content. A previous study on natural perennial ryegrass populations also revealed that even though natural populations showed on average lower agronomic performances, their range of variation was much larger (Balfourier and Charmet 1991). Such results indicate that natural diversity holds potential for use in breeding programs to further improve the agronomic value of cultivars as suggested by other authors (Charmet et al. 1989). Novel combinations of natural sources could likely generate original phenotypes including more extreme trait values and/or new combinations of trait values that would be worthwhile for both climatic adaptation and agronomical value.

Valorising GWAS results

Population allele frequencies and mean phenotypes per population have already been successfully combined for association studies to detect major effect loci in various species such as Drosophila melanogaster (Bastide et al. 2013), humans (Gaj et al. 2012; Riaz et al. 2016) and perennial ryegrass (Fè et al. 2015a). This was also the case in the present study in which areas of the genome including causal genes for various traits were revealed (Table S9).

A notable example is a region on chromosome 7 which appeared linked to heading date as well as spring growth, leaf lamina width, spring forage quality traits, flowering stem density and straw height (Figure 3). Another region on chromosome 4 appeared linked to first year heading, heading date and aftermath heading in multiple environments. One hypothesis is that such regions of the genome contain pleiotropic genes (Stearns 2010) meaning that all associated traits would have a common genetic determinism. Another hypothesis is that such regions contain multiple genes in linkage (Zeven and Harten 1979) that would correspond to co-adapted functions and for which limited allelic arrangements provide high adaptive values.

The number of QTL detected by GWAS and the percentage of phenotypic variance explained varied greatly between traits (Figure 4 and Table S9). The effect of a single SNP on phenotypic variation in a GWAS model which accounted for kinship was likely underestimated due to the strong confusion between the genetic background and the distribution of values for most traits. As such the effect of each SNP detected by GWAS (which accounted for kinship) was also computed using a linear regression model between trait values and SNP allele frequencies which did not account for kinship. That method on the contrary likely overestimated the SNP effect. The common GWAS issue of missing heritability (Manolio et al. 2009) was more or less important depending on the trait. In no case could the SNP markers explain all of the inter-population genetic variance.

The four following reasons could explain the poor capacity of GWAS to detect significant makers for a number of traits. (i) Even if the number of markers used was high, some QTL could have been missed due to the fast LD decay. (ii) Because of confusion between genetic kinship and phenotypic similarity, accounting for kinship in GWAS models could have hampered the detection of QTL. For natural diversity in which diversifying selection and gene flow are strong, it has been proposed that only a small number of large effect QTL that shape inter-population trait variability may show high levels of allelic differentiation (Corre and Kremer 2003; Storz 2005; Bell 2010). A large portion of the loci involved in the variability of a trait may behave like neutral markers in terms of allelic differentiation. Consequently, GWAS models accounting for kinship could only detect large effect adaptive QTL. Furthermore, the continuous genetic structure due to isolation by distance has for consequence that few or no significant SNPs could be detected by GWAS for traits showing high spatial autocorrelation. (iii) A high proportion of the GBS SNP markers had a rare minor allele (present in few populations and/or at low frequency in populations). As reported from previous studies (Sanchez-Bermejo et al. 2015), rare alleles may have been responsible for a substantial part of the phenotypic variation of many traits and make their genetic determinism difficult to apprehend. (iv) Another hypothesis is that the variation of some traits results from many genes with small effects. It was suggested (Nicolae et al. 2010; Zhao et al. 2011) that the integration of transcriptome data could improve the ability to detect moderate effect loci if gene expression levels can be associated to phenotypic variation, i.e., in the case of expression quantitative trait loci (eQTL).

For all traits, there are cases of phenotypically similar populations that are spatially and genetically distantly related even if genetic distance is calculated using the allele frequencies of the SNPs most significantly associated to the trait (Figure S4). This suggests that the causal polymorphism can vary according to geographical regions. This shows that the search for natural genetic resources to improve a given trait, or adaptation to a given climatic constraint, should generally not focus on a single restricted geographical region. Original genotypes recombining valuable alleles at various loci could be designed through transgressive segregation (Rieseberg et al. 1999).

Population × environment (G×E) interactions should be taken into account when implementing GWAS. Many phenotype-SNP associations were environment-dependent (Table S10), a likely common issue in association studies that is often ignored (Hirschhorn et al. 2002). The effect of environment on biological systems is multifactorial and it is thus difficult to determine which environmental factor, or combination of factors, interacts with the QTL. Therefore, it may be difficult to predict the effect of environment interactive SNPs in a new environment and so to use such markers reliably. However, when a SNP marker is located within a gene whose function is known, assumptions may be formulated to explain SNP × environment interactions. For example in this study, the effect of a SNP in the LpVRN2 vernalisation gene (Andersen et al. 2006) varied according to year and location (Figure 6) as if cold winters were less likely to lead to differences in heading date determined by variation in vernalization response. To confirm such relations and to further investigate GxE interactive traits, it may be necessary to perform association studies using the reaction norm slope of a given trait to an environmental variable obtained across many different environments (Sanchez-Bermejo et al. 2015) .

Genomic predictions

The results of this study are very promising as to the ability of GP models to accurately predict the value of natural populations of a grassland species for many traits relevant to agronomical value and/or environmental adaptation (Figure 7, Table S7). GP were already tested on various forage species notably alfalfa for which Jia et al. (2018) found the highest accuracy (0.65) for autumn plant height, switchgrass for which Ramstein et al. (2016) did not find accuracies exceeding 0.5 and perennial ryegrass for which results are very variable as described in the following. Faville et al. (2018) reported a study in which 517 mother plants were sampled from a commercial breeding program, genotyped for over one million SNP markers and their half-sib progeny families were phenotyped. That study led to predictive abilities that did not exceed 0.52 for heading date. Pembleton et al. (2018) reported a study based on 18 years of phenotypic data recorded in a commercial breeding program involving 714 synthetic populations that were genotyped for over 28,000 SNP markers. They reported a predictive accuracy of 0.31 for biomass yield and 0.76 for heading date. Other than the relatively high number of available SNP markers and phenotyped accessions, what distinguishes our study is both the use of allele frequencies of populations rather than individual genotypes and the use of natural highly diverse genetic material.

More or less ancient selection may have been the main driving force that shaped the extant phenotypic variability of natural populations (Rieseberg et al. 2002). The natural populations we studied originated from contrasting environments likely requiring different local adaptive trait optima and imposing different selection pressures at the underlying QTL. These QTL should thus have outstanding allelic distributions as compared to neutral loci affected mainly by drift and migration. Yet, due to the outbreeding nature of perennial ryegrass, gene flow between spatially close populations is likely fairly strong, as suggested by the correlation of 0.44 between geographic and genetic distances which suggests isolation by distance. As such, it is likely that high levels of allele differentiation from the genome-wide genetic structure only exist for the few largest effect QTL under strong selection (Corre and Kremer 2003; Bell 2010) . The allele distributions of most loci involved in the phenotypic variability of adaptive traits are likely confounded with the neutral continuous genetic structure. For non-adaptive traits, the distribution of population values and of allele frequencies at most of the underlying QTL is expected to follow neutral patterns and to be more specifically ruled by isolation by distance. As such, for most traits the patterns of isolation by distance for the phenotypic values and for the genetic background are comparable, thus making kinship analyses (genetic relatedness) a very powerful tool for predicting phenotypic differences at the ‘species scale’. The high predictive abilities obtained for phenotypic traits could be largely thanks to a good estimation of kinship with high density molecular markers.

The preceding hypothesis could explain the differences of predictive ability between natural populations and cultivars (Table S7), which could be related to the genetic architecture of predicted traits. For example, heading date is thought to be under the control of few large effect genes (Fè et al. 2015a) as our GWAS results confirmed, whereas autumn growth appeared more likely to be under the control of many small effect genes. It could thus be that heading date was well predicted both for natural populations and cultivars because the prediction model uses information linked to loci that have a large and stable effect regardless of the genetic background. Conversely, it could be that autumn growth, which shows strong spatial autocorrelation, was well predicted for natural populations largely because of a good estimation of kinship but the models were far less effective for the cultivars for which the genetic structure has been broken by genetic recombination cycles.

The use of GP models to describe genebank accessions was already evaluated for various domesticated plant species such as wheat (Crossa et al. 2016), sorghum (Yu et al. 2016), cauliflower (Thorwarth et al. 2018), pea (Burstin et al. 2015), soybean (Peixoto et al. 2017) and for the natural diversity of Arabidopsis (Kooke et al. 2016). In most of these studies, as well as in ours, predictive abilities were remarkably high for many traits and genetic structure (continuous or stratified) seemed to play an important role. For example, Yu et al. (2016) reported an accuracy of 0.84 when predicting dry biomass yield for 200 accessions of sorghum with a model that used over 250,000 SNP markers. However, the high accuracy appeared to be mainly due to the ability of the model to reveal phenotypic differences between genetic groups and not so much within groups. Population structure is expected to be high in genebank collections capturing large portions of intra-specific variability compared to that within commercial breeding material (Crossa et al. 2016). Such structure could partly explain the high predictive abilities observed in the above quoted studies in which models likely used kinship to predict breeding values of traits whose distribution was more or less confounded with genetic structure. Here, we argued that the continuous genetic structure in our genetic material, as previously reported by Blanco-Pastor et al. (2019), likely explains the high predictive abilities we obtained.

Practical use of GP for mining natural genetic resources

Describing a non-phenotyped portion of genetic resources by the mean of GP models calibrated with phenotypic data recorded on a another portion of these genetic resources leads to query about which resources should be phenotyped and what genotyping effort should be made.

Optimization of calibration sets for GP models, based on stratified sampling was previously shown to be effective (Isidro et al. 2015). The methods used in the present study were similar except that the size of the desired calibration set equals the number of clusters to choose from. If all populations or individuals from a genetic resource collection are genotyped, our results show that genetic distance optimization can be successfully implemented. However, if genotyping information is not available, our results show that it can be efficiently substituted by geographical information if phenotypic diversity is spatially structured. Stratified sampling in spatial clusters was already shown to be an efficient way to choose a core collection of natural populations of perennial ryegrass (Balfourier et al. 1998). The best choice between optimization based on genetic or geographical information depends on the trait and the relative size of the calibration set. The use of stratified sampling after clustering based on both genetic similarity and spatial proximity appeared to be the most stable (across traits and calibration set sizes) procedure to choose a calibration set. Such a procedure makes possible the choice of a single subset of accessions that can be used to calibrate prediction models for any trait without prior knowledge of their heritability. This is a benefit compared to method such as CDmean optimization (Rincent et al. 2012) which gave us similar results in terms of predictive ability gains relative to random sampling (not shown) but must be applied trait by trait.

Regarding the genotyping effort, it could be considered to optimize the number and choice of SNP markers to use. Many previous studies showed, as we did, that the predictive ability of genomic models improves as the number of included SNP markers increases until a saturation is reached (Spindel et al. 2016; Thorwarth et al. 2018). For any trait, 1000 random SNPs appears sufficient to reach near maximal predictive abilities, which suggests that the estimation of kinship does not indefinitely improve with the addition of extra genetic information. Furthermore, the GP of many traits was improved by including markers identified as most significant in a GWAS based on the calibration set data. Such markers may not necessarily be linked to QTL but may simply have an allele whose population frequency distribution is particularly well correlated to the phenotypic diversity for the trait considered.

In the present study, the populations were sampled across the natural geographic range of perennial ryegrass so that most of the genetic diversity of the species may be represented. Hence, if GP models calibrated on a sub-set of the sampled populations can accurately predict the remaining populations, they should also accurately predict any other natural genetic material already present in a genebank or collected in its native environment. Consequently, this study shows that GP can greatly improve capabilities to mine the whole natural genetic diversity of a species for many traits of agronomical and/or ecological importance. Furthermore, the objective of choosing an optimized set of populations for the calibration of GP models, as presented in this study, is quite similar to the objective of choosing an optimized core collection representative of the genetic diversity to be safeguarded in a genebank. Setting up such core collections is relevant to maximize the ratio of preserved diversity to conservation costs but also to study the diversity of a species with a representative sample of limited size (van Hintum et al. 2000).