The estimation of genomic heritability was based on the following bivariate statistical model across environments:y=[y1y2]=Xβ+Zu+e,(1)where y1, y2 are the phenotypic values (absorbance for a given wavelength) in each environment, β is a vector of fixed environment effect, u is a vector of random polygenic effect with var(u)=[σu12σu12σu21σu22]⊗K, K being the scaled realized relationship matrix (see below), e is a vector of independent and normally distributed residuals with var(e)=[σe1200σe22]⊗I, X and Z are design matrices relating observations to the effects. SNPs were used to estimate the genomic relationship matrix (A) between individuals, following the formula of VanRaden (VanRaden 2008)A=1L∑l=1L(Gi,l−pl)(Gj,l−pl)σ2,(2)where Gi,l and Gj,l are the genotypes of individuals i and j at marker l (G.,l=0 or 1 for homozygotes, 0.5 for heterozygotes), pl is the frequency of the allele coded 1 for the marker l, and σ2 is the average empirical marker genotype variance. K was obtained by scaling A to have a sample variance of 1 (Kang et al. 2010; Forni et al. 2011). Genomic heritability was estimated for each wavelength within each environment (m) as follows:hm2=σ^um2σ^um2+σ^em2,(3)with σ^um2 and σ^em2 the REML estimates of σum2 and σem2, obtained with the Newton-Raphson algorithm implemented in the R package sommer (Covarrubias-Pazaran 2016). Following Yamada et al. (1988), the variance/covariance estimates from the previously defined bivariate mixed-model were used to compute estimates of genetic (σ^G2), genetic by environment (σ^G×E2) and residual (σ^ϵ2) variances across sites as follows: σ^G2=σ^u12, σ^G×E2=12(σ^u12+σ^u22)−σ^u12, and σ^ϵ2=12(σ^e12+σ^e22).

Association mapping was carried out along spectra considering the absorbance at a given wavelength as a bivariate trait (corresponding to the two environments) and using previous estimates of genetic and residual variances (EMMAX philosophy as previously proposed in the multi-trait mixed-model approach (Korte et al. 2012)).

The efficiency of genomic and phenomic predictions was evaluated by cross-validations in two types of scenarios (Figure 1). In scenario S1, NIRS analysis and cross-validation were applied to the same environment (Figure 1 a). In scenario S2, cross-validation was applied to independent environments: the environment(s) in which NIRS was collected and the environment in which the cross-validation was applied (calibration and prediction) were different (Figure 1 b). In S1, the objective was to limit expensive or labor-demanding phenotyping to a calibration set of reduced size and to predict the remaining individuals using NIRS. In scenario S2, one experiment (or a nursery) was dedicated to collecting the NIRS of the calibration set and the predicted set, and a multi-environment trial was dedicated to phenotyping the calibration set. The main difference between S1 and S2 was that in S1, we can expect NIRS to be more related to the phenotypic data than in S2.

For both scenarios, 5- and 8-fold cross-validation procedures repeated 20 times were used for poplar and wheat, respectively. A larger fold-number was considered for wheat in comparison to poplar because the sample size in the wheat dataset (n = 223 in the panel, and n = 161 in the subpanel) was lower than the sample size in the poplar dataset (n = 562). We would like to emphasize that in both scenarios, the predicted set was only characterized by genotypic and NIRS (obtained in the reference site) data. In particular in scenario S2, the only information available from the predicted environment is the phenotypic data of the calibration set (Figure S1). Predictive ability was computed as the Pearson correlation between the predictions and adjusted means. For genomic predictions, we tested two complementary reference models: G-BLUP and Bayesian LASSO (Park and Casella 2008; de los Campos et al. 2009). The underlying assumptions of these two models are that the SNP effects are normally distributed for G-BLUP, whereas Bayesian LASSO allows for departure from normality (i.e., SNPs with bigger effects). G-BLUP and Bayesian LASSO were run with the R packages rrBLUP (Endelman 2011) and BGLR (de los Campos et al. 2013), respectively. For Bayesian LASSO, the chain was composed of 30,000 iterations with a burn-in of 5,000 iterations, and the hyperparameter λ was chosen as recommended in de los Campos et al. (2013). For phenomic predictions, we used RR-BLUP but considered NIRS data instead of molecular markers. Prior to the analysis, the pretreated NIRS matrices were centered and scaled for each wavelength. For wheat, we also tested to include the spectrum of the two tissues (or of the same tissue but collected in the two environments) in the statistical model. For this we simply joined the two matrices of spectrum into one single matrix.

We ran simulations to illustrate the expected genetic gain with GS and PS that would be achieved in one cycle of selection for various combinations of costs and reliabilities. Reliability was defined as the squared correlation between true breeding values (TBV) and the genomic or NIRS predicted values (PV).

We considered a situation in which a given budget (200,000 €) was available to predict the performances of selection candidates with NIRS or genotyping. Depending on the costs of the methods (DNA extraction and genotyping for GS or tissue sampling and NIRS acquisition for PS), we computed the number of selection candidates (N) that could be analyzed. The TBV and genomic or NIRS PV of these N individuals were then sampled from a multivariate normal distribution with means equal to 0, variances equal to 1 and covariance equal to the square root of reliability (R package mvtnorm (Genz et al. 2018)). The expected genetic gain was then computed as the difference between the average TBV of the 400 individuals having the best PV and the average TBV of the population (equal to 0). We selected 400 individuals because for many species, it is feasible to apply heavier phenotyping (multi-environment trials) on a few hundred individuals. We considered two situations; in the first situation, the expected genetic gain of GS and PS was computed for various genotyping and NIRS costs with a reliability set to 0.4. In the second situation, the reliability of GS and PS varied between 0.3 and 0.6, and genotyping and NIRS costs were set to 50 € and 4 €, respectively. For each combination of parameters (reliabilities and costs of GS and PS), the simulation procedure was repeated 1000 times to obtain stable results. Because genotyping and NIRS costs are highly dependent on the species and the number of samples analyzed, we let the genotyping costs (DNA extraction and genotyping itself) vary between 25 € and 100 € and the NIRS costs (sample treatment and NIRS analysis itself) vary between 1 € and 8 € in the first situation.

To provide concrete examples, we applied this simulation process with the reliabilities and costs that we experienced for wheat and poplar. GS costs were between 35 € and 50 € per individual for wheat and poplar, respectively, and PS costs were between 3 € and 2.5 € per individual for wheat and poplar, respectively. Reliabilities were estimated as the square of predictive abilities estimated by cross-validation divided by the heritability of the adjusted means. For each combination of trait, scenario, and NIRS data considered (tissue, environment), the increase in expected genetic gain using PS instead of GS was computed with the best performing GS model as a reference.

The datasets generated during and/or analyzed during the current study are available in the INRA Dataverse repository (https://data.inra.fr/). They can be accessed with the following link http://dx.doi.org/10.15454/MB4G3T. R functions used for comparing the predictive ability of SNP and NIRS through cross-validations have been deposited on github (https://github.com/visegura/PS). Supplemental material available at Figshare: https://doi.org/10.25387/g3.7243256.

We first sought to characterize the ability of NIRS to capture genetic variability by estimating genomic heritability and partitioning the variance into genetic (G), genetic by environment (G×E) and residual variances (ϵ) along the NIR spectrum collected on a panel of winter wheat (leaves and grains) and a population of black poplar (wood) grown in two contrasting environments. For both species and tissues, genomic heritability was highly variable along the spectrum with peaks above 60%, showing the existence of strong polygenic signals for some wavelengths (Figure 2, Figure S2, and File S1). For a given species, the proportion of G×E variance was on average across all the wavelengths equal to 8% (poplar), 10% (wheat leaves) or 16% (wheat grains) (Figure 2 and File S1). It is interesting to note that for at least half of the wavelengths, the cumulative proportion of G and G×E variances was above 15%, showing that the NIR signal was often partially related to genetics. The kind of tissue analyzed by NIRS seemed to matter, as shown by the comparison of variance partition along spectra obtained on wheat leaves and grains. G and G×E variances were higher and more stable along the spectrum for grains than for leaves.

We ran association mapping along the NIR spectrum to identify wavelengths associated with major QTL (Figure S2). In poplar, the signal appeared to be mainly polygenic with very few QTL detected, and the largest SNP R2 was below 0.025 for any wavelength. In contrast, in winter wheat, we detected numerous large-effect QTL. For some wavelengths, a single SNP could have an R2 of 0.23 for leaves and of 0.11 for grains, and this SNP could be in spectrum regions of high or of low genomic heritability (Figure 2). This finding means that depending on the wavelength, NIRS could capture highly polygenic relationships (wavelengths with high genomic heritability) or could tag specific regions of the genome (major QTL). These two kinds of wavelengths can be useful for making predictions because they can potentially track the two main factors responsible for GS accuracy: relatedness and linkage disequilibrium.

We estimated the efficiency of GS and PS to predict the performance of new individuals within a cross-validation framework. The performances of the individuals in the validation set were predicted with genotypic information in GS (G-BLUP and Bayesian LASSO models) and with NIRS only in PS (RR-BLUP model). We considered two scenarios: in S1, NIRS analysis and cross-validation were performed in the same environment (Figure 1 a), whereas in S2, the environments in which the cross-validation was applied were different from those in which NIRS was obtained (Figure 1 b). The broad-sense heritabilities of the adjusted means were above 0.8 for all traits in each environment (Table S1).

In wheat, the predictive abilities of PS were highly variable and appeared to be dependent on the predicted trait and on the environment and tissue in which NIRS was measured (Figure 3 a, b, c, d and Figure 4). While combining NIRS collected in different environments or different tissues increased the predictive ability, this increase did not occur systematically. One major result is that for both traits, NIRS could lead to better predictions than molecular markers, even in the six independent environments (Figure 3 c, d and Figure 4). The gain with NIRS in comparison to molecular markers in S2 was up to 34% and 22% for heading date and grain yield, respectively. In each S2 environment and for both traits, there was always a type of NIRS that performed better or as well as the best GS model (Figure 4). The gain was even stronger in S1: NIRS led to an increase in predictive ability of up to 53% and 117% for heading date and grain yield, respectively. In poplar, the predictive abilities with NIRS were always lower than those with SNP, except for growth traits under S1 (Figure 3 e). In the other cases, the predictive ability with NIRS varied depending on the trait and scenario considered, but they were always significantly greater than 0. In general, they were higher when the spectra were collected in the same environment (S1) than when spectra from another environment were used (S2), except for bud flush evaluated in one site and bud set evaluated in another site. Interestingly, irrespectively of the scenario, for some traits apparently unrelated to wood chemical properties, such as resistance to rust or bud set, NIRS predictive abilities were fairly high ranging between 0.34 and 0.53.

To further evaluate the potential of PS with respect to GS, the expected genetic gain with both approaches was compared in a simple scenario in which a budget of 200,000 € could be spent to genotype or analyze the NIRS of selection candidates. The difference in efficiency between GS and PS was highly dependent on the genotyping and NIRS costs and on the reliability of the two approaches (Figure 5). In the scenarios that we considered here, the expected gain of using PS instead of GS was between 11% and 127%. In extreme scenarios in which genotyping was cheap (25 €) and NIRS was expensive (8 €) or in which GS reliability (0.6) was much higher than PS reliability (0.3), PS was still better than GS. We applied the simulation process with the reliabilities and costs obtained in the wheat example (35 € for genotyping and DNA extraction and 3 € for sample treatment and NIRS acquisition). The increase of expected genetic gain with PS in comparison to GS was between +60% and +127% for heading date and between -10% and +222% for grain yield, depending on the tissue and environment used for NIRS acquisition and scenario considered (Table S2). In poplar, considering genotyping and NIRS acquisition costs of 50 € and 2.5 €, respectively, as well as the reliabilities estimated with cross-validation predictive abilities, the expected gain in genetic progress varied depending on the trait and scenario considered (Table S3). It was mainly positive for growth traits (-2–93%), bud set (-6–25%) and rust resistance (-10–21%), whereas for bud flush, NIRS prediction did not seem to provide any advantage over regular SNP-based prediction.