The dataset included a total of 3,771 bread wheat lines evaluated at the Campo Experimental Norman E. Borlaug in Ciudad Obregón, México over the course of four breeding cycles: 2013-14, 2014-15, 2015-16, and 2016-17. In each cycle, lines were sown under five differentially managed treatments: Optimal Bed, Optimal Flat, Moderate Drought, Severe Drought, and Heat. Descriptions of the managed treatments are given in Table 1. The 20 managed treatment-breeding cycle combinations are herein referred to as site-years. Within each of the five managed treatments, 1,092 lines were arranged into 39 trials in a α-lattice design with three replicates and six incomplete blocks per replicate. Each replicate contained two repeated checks. ‘Kachu #1’ was sown as a check in the Optimal Bed, Optimal Flat, and Moderate Drought managed treatments. ‘Baj #1’ was used as a check in the Severe Drought and Heat managed treatments. ‘Borlaug100 F2014’ was sown as the other check in all managed treatments. All lines within a breeding cycle were evaluated in all five managed treatments, while no lines other than checks overlapped between cycles. Records for some lines were removed from the analysis due to unavailability of genotypic or agronomic data, resulting in the final dataset: 588 lines in 2013-14, 1,033 lines in 2014-15, 1,063 lines in 2015-16, and 1,087 lines in 2016-17, for a total of 56,565 plots phenotyped. Each breeding cycle contained full-sib families with an average of two full-sibs per family.

Grain yield (GY) in t ha^-1 and lodging evaluated on an ordinal scale (0: no lodging; 5: completely lodged) were assessed in all three replicates. Heading date was recorded as the date of 50% spike emergence within the plot. Maturity date was recorded as the date of senescence of the peduncle for 50% of stems in the plot. Days to heading (DTHD) and days to maturity (DTMT) were measured as the number of days to reach heading and maturity, starting from the date of the first irrigation if sowing was on dry soil or from the sowing date if sowing was in pre-irrigated fields. DTHD and DTMT were assessed for the first replicate only due to the high heritability of those traits.

Hyperspectral reflectance data were collected with a hyperspectral camera (A-series, Micro-Hyperspec VNIR, Headwall Photonics, Fitchburg, MA, USA) as part of the Alava Remote Sensing Spectral Solution (ARS3, Alava Ingenieros, Madrid, Spain) mounted in a Piper PA-16 Clipper aircraft. The camera’s sensor had a 12-bit radiometric resolution, covering the light spectrum in the 380–850nm region with a 7.5nm full width at half maximum and set with an integration time of 18ms. Spectral binning resulted in 62 wavelengths between 398-847nm (Rodrigues et al. 2018). The flights were scheduled around noon (GMT-7) and aligned to the solar azimuth angle at a height of 300m, resulting in 30cm resolution. The aerial imagery acquisition during the growing season was spaced approximately at seven- to ten-day intervals on mostly clear days.

Radiometric calibration of the sensor was done using coefficients derived from a calibrated uniform light source and an integrating sphere (CSTMUSS2000C Uniform Source System, LabSphere, North Sutton, NH, USA). Dark frame correction was performed for each flight dataset. Atmospheric calibration was performed using irradiance measurements acquired at the beginning and end of each flight using a Jaz spectrometer with a CC-3 Cosine Corrector (Ocean Optics Inc, FL, USA) for the 2016-17 breeding cycle. For remaining three breeding cycles, irradiance was modeled using aerosol optical depth from sun-photometer measurements (Microtops II, Solar Light Company, Glenside, PA, USA) based on the SMARTS simulation model (Gueymard 1995; Gueymard 2005).

Ortho-rectification and georeferencing of the imagery were performed using PARGE (ReSe Applications Schläpfer, Wil, Switzerland) based on data from the inertial navigation system (INS) attached to the camera (IG-500N model, SGB systems S. A. S., Carrières-sur-Seine, France). Hyperspectral reflectance data were extracted from the aerial images using the mean value of the pixels inside the central area of each observed plot, avoiding 0.5m from the plot border.

Each hyperspectral phenotyping time-point within each site-year was assigned a developmental growth stage classification according to the predominant growth stage of the lines within the site-year at the time of phenotyping (Table 2). The vegetative (VEG) stage was defined as the period between germination and 50% of plots at heading. The heading stage (HEAD) was defined as the period between 50% of plots at heading and 100% of plots at heading. The grain fill stage (GF) was defined as the period between 100% of plots at heading and 100% of plots at maturity.

All lines were genotyped using genotyping-by-sequencing (Elshire et al. 2011) according to the pipeline described in Poland et al. (2012). From the initial set of 34,900 single nucleotide polymorphisms (SNPs), 8,519 SNPs remained after excluding all markers with more than 70% of missing data or minor allele frequency less than 0.05. For each marker, missing data were imputed using the sample mean of observed values (Poland et al. 2012).

Within each site-year, best linear unbiased estimates (BLUEs) were calculated for the agronomic traits and for each hyperspectral band at each phenotyping time-point using the following model:

where y_ijkl is the trait value for genotype i within trial j, replicate k, and block l; μ is the overall mean; g_i is the fixed effect for genotype i; t_j is the random effect for trial j, which are assumed to be independently and identically distributed according to a normal distribution with mean zero and variance σ_t², that is, t_j ∼iid N(0, σ_t²); r_k(j) ∼iid N(0, σ_r²) is the random effect for replicate k within trial j; b_l(jk) ∼iid N(0, σ_b²) is the random effect for block l within replicate k and trial j; and ε_ijkl ∼iid N(0, σε2) is the residual effect. For DTHD and DTMT, which were recorded for only one replicate, the replicate and block effects were excluded from the model. For validation of prediction models, best linear unbiased predictions (BLUPs) for GY within each site-year were calculated by treating the effect for genotype g_i in model (1) as random with g_i ∼iid N(0, σ_g²) and by including a covariate for lodging as a fixed effect, fit only when lodging was observed within a site-year. To correct for the influence of phenology, BLUPs for GY were calculated again using DTHD included as a fixed effect in model (1).

BLUEs were also calculated for each hyperspectral band on a developmental growth stage basis according to the growth stage classifications detailed in Table 2. BLUEs for each hyperspectral band for the VEG, HEAD, and GF stages were estimated by fitting the following model:

where y_ijkl is the hyperspectral reflectance trait value for genotype i on time-point j and within trial k, replicate l, and block m; μ is the overall mean; g_i is the fixed genetic effect for genotype i; d_j ∼iid N(0, σ_d²) is a random effect for time-point j; t_k(j) ∼iid N(0, σ_t²) is random effect for trial k nested within time-point j; r_l(jk) ∼iid N(0, σ_r²) is the random effect for replicate l nested within time-point j and trial k; b_m(jkl) ∼iid N(0, σ_b²) is the random effect for block m nested within time-point j, trial k, and rep l; and ε_ijklm ∼iid N(0, σε2) is the residual effect. BLUEs were also calculated for each hyperspectral band across all available phenotyping time-points (ALL) using model (2).

Broad-sense heritability within site-years was calculated for GY (Bernardo 2010) and each hyperspectral band as:

where σ_g² is the genetic variance, σε2 is the error variance, and nreps is the number of replicates (nreps = 3). σ_g² and σε2 were derived by fitting model (1) with the effect for genotype g_i treated as random with g_i ∼iid N(0, σ_g²). For each breeding cycle, Pearson’s correlations for GY between managed treatments were calculated using the BLUEs derived from model (1). Pearson’s correlations were also calculated between GY BLUEs and the BLUEs for each hyperspectral band calculated in model (2).

Genetic relationships between individuals were modeled using genomic markers, pedigrees, and hyperspectral reflectance phenotypes. The genomic relationship matrix (G) was calculated according to Endelman and Jannink (2012). The additive relationship matrix (A) was derived from pedigrees and calculated as twice the coefficient of parentage.

Hyperspectral reflectance-based relationship matrices (H) were calculated within each site-year using the hyperspectral BLUEs calculated for 1) the individual time-points from model (1), 2) the developmental growth stages from model (2), and 3) all time-points from model (2). The matrices were calculated as H=SS′nbands, where S is a matrix of the centered and standardized BLUEs of the hyperspectral bands and nbands is the total number of hyperspectral bands observed (nbands = 62).

Processing of the hyperspectral images was performed using the ARS3 hyproQ software (Álava Ingenieros, Madrid, Spain). Plot polygons for tabular data extraction were generated using ArcGIS (ESRI, Redlands, California, USA). Images were aligned manually in ArcGIS if they did not overlay the plot polygons due to INS inaccuracy.

All models were fit using the R statistical programming language (R Core Team 2018). Basic models (1-3) were fit with the package “ASReml-R” (Gilmour et al. 2014) for R, while the prediction models (4-7) were fit using the “BGLR” package (de los Campos and Pérez-Rodríguez 2014) for R. The marker-based genetic relationship matrices were calculated using the A.mat() function within the “rrBLUP” package (Endelman 2011) for R. The coefficients of parentage for the pedigree relationship matrices were estimated using the “Browse” application within the International Crop Information System software package (McLaren et al. 2000).

The above models were fit to assess model prediction accuracy for three prediction strategies representing different testing and evaluation problems relevant to plant breeding programs: 1) within site-year, 2) within breeding cycle/across managed treatments, and 3) across breeding cycles/within managed treatment.

A train-test (TRN-TST) validation scheme (Daetwyler et al. 2013) of 20 random partitions was used to assess model prediction accuracy for all prediction strategies. For within site-year prediction, a random 80% of records within a given site-year were assigned to the TRN set, and the remaining 20% were used as the TST set for prediction. This reflects a scenario in which prediction of the trait values for a set of unobserved lines is performed within a managed treatment and breeding cycle of interest. Models (4) and (5) were fitted without the E_i term for site-year and models (6) and (7) were not tested for the within site-year prediction scheme because multiple site-years were not considered.

For prediction within breeding cycle/across managed treatments, predictions were carried out across the five managed treatments within a single breeding cycle. The TRN set consisted of all records from four out of the five managed treatments within a breeding cycle plus 20% of records from the fifth managed treatment. The remaining 80% of records from the fifth managed treatment were assigned to the TST set for prediction. In this situation, lines have been previously characterized in all managed treatments except one, where only 20% of records were available.

For prediction across breeding cycles/within managed treatment, predictions were performed across the four breeding cycles but within a single managed treatment. The TRN set contained all records from three out of the four breeding cycles for a given managed treatment plus 20% of records from the fourth breeding cycle. The remaining 80% of records from the fourth breeding cycle were assigned to the TST set for prediction. This scenario predicts performance of unobserved lines that are 1-3 breeding cycles removed from lines in the training set.

Table S1 illustrates examples of TRN-TST partitioning for the three prediction schemes. For all prediction strategies, model prediction accuracy was assessed as the Pearson’s correlation between predictions for GY and GY BLUPs with and without correction for DTHD calculated in model (1). Reported values are the mean and standard deviation of the 20 random TRN-TST partitions implemented. The same partitioning “splits” were used to assess all models to ensure fair comparisons.

All phenotypic and genotypic data required to confirm the results presented in this study are available on CIMMYT Dataverse link: hdl:11529/10548109. The “GID” column denotes the unique identifiers for the genotypes. Supplemental material available at Figshare: https://doi.org/10.25387/g3.7653473.

To evaluate the potential of integrating aerial hyperspectral reflectance phenotypes into GS to improve prediction accuracy for GY in wheat, we deployed an airplane equipped with a hyperspectral camera to phenotype wheat breeding trials in 20 site-years at multiple time-points throughout the growing season. Across the site-years, the Optimal Bed and Optimal Flat managed treatments had higher GY than the stressed Moderate Drought, Severe Drought, and Heat managed treatments (Figure 1). The standard deviations for GY remained relatively stable across site-years, ranging from 0.30 to 0.68 t ha^-1 (Table S2). Broad-sense heritability for GY was high across site-years, ranging from 0.58 to 0.94 (Table S2). Overall, correlations for GY between managed treatments were moderately positive (Figure S1). The managed treatments that received similar levels of irrigation (e.g., Optimal Bed with Optimal Flat; Moderate Drought with Severe Drought) showed significant correlations (p-value ≤ 0.05) in all breeding cycles.

Between 3 and 11 hyperspectral phenotyping time-points were collected within each site-year (Table 2). In most site-years, at least one hyperspectral phenotyping time-point was collected during each of the three developmental growth stages. While the managed treatments within a breeding cycle were often phenotyped on the same date, the growth stage at the time of phenotyping frequently differed among the managed treatments (Figure 2). During the 2015-16 cycle, technological issues with the camera prevented early season phenotyping. As a result, no hyperspectral data were collected during the VEG stage. The Optimal Bed, Optimal Flat, and Heat managed treatments had hyperspectral phenotypes for the HEAD and GF stages only, while the Moderate Drought and Severe Drought managed treatments had hyperspectral phenotypes for the GF stage only due to early maturation of the crop under water stress.

Broad-sense heritability estimates for the hyperspectral bands were generally between 0.5 and 0.8 for most phenotyping time-points (Figure 3). Heritabilities were relatively homogeneous across individual time-points and developmental growth stages. In all site-years, the lowest heritabilities were observed for the hyperspectral bands between 398 and 425 nm. Heritabilities in the 2013-14 Optimal Flat site-year were notably lower than in other site-years, falling between 0.00 and 0.25, though the GF time-points were more heritable with values around 0.40. Low heritabilities (<0.20) were also observed for the 9 March time-point in the 2015-16 Moderate Drought site-year.

Correlations between the individual hyperspectral bands and GY ranged from -0.68 to 0.64 across the 20 site-years and were most frequently the strongest above 500 nm (Figure 4). While there were some observable similarities in correlation patterns among time-points taken during the same growth stage within the same site-year, patterns were not uniform across site-years, breeding cycles, or managed treatments. In some site-years, there were clear differences in correlation patterns between developmental growth stages. For example, in the 2014-15 Optimal Bed site-year, correlations between GY and hyperspectral reflectance in the 398-700 nm range were around 0 during the VEG stage but became progressively more negative during the HEAD and GF stages, reaching around -0.4. For the same managed treatment during the 2016-17 breeding cycle, distinct differences in correlation patterns between the three developmental growth stages were also observed; however, the patterns do not reflect those observed in 2014-15. Despite irregularities in correlation patterns, the correlations between hyperspectral bands from the GF stages and GY were, on average, stronger by 0.10 than for hyperspectral phenotypes taken during the VEG and HEAD stages.

In this study, we have proposed a multi-kernel GBLUP model that uses genomic marker-, pedigree-, and hyperspectral reflectance-derived relationships matrices to model the genetic main effects and G×E interactions across environments within a bread wheat breeding program. We have shown that deriving relationship matrices from aerial hyperspectral reflectance phenotypes can effectively predict GY in wheat within and across managed treatments and breeding cycles. Accuracies when testing single-kernel models using hyperspectral reflectance data alone are similar to those achieved with markers or pedigrees. Our results also show that combining markers/pedigrees with hyperspectral reflectance data in multi-kernel models can increase accuracies over single-kernel approaches, but in some prediction scenarios, these increases were modest. We also suggested a method for addressing the issue of unbalanced HTP datasets involving the classification of time-points according to the predominant developmental growth stage observed at the time of phenotyping. Accuracies of multi-kernel models were roughly equivalent irrespective of the growth stage in which hyperspectral phenotyping was performed. Further research on how best to leverage multi-temporal phenotypes when the amount of data differs across site-years is needed. The methods we have proposed provide a simple and computationally efficient approach for integrating highly dimensional aerial HTP information into genomic selection and should be tested on other forms of high dimensional biological data.

All phenotypic and genotypic data required to confirm the results presented in this study are available on CIMMYT Dataverse link: hdl:11529/10548109. The “GID” column denotes the unique identifiers for the genotypes. Supplemental material available at Figshare: https://doi.org/10.25387/g3.7653473.