Aerial High-Throughput Phenotyping Enabling Indirect Selection for Grain Yield at the Early-generation Seed-limited Stages in Breeding Programs

Breeding programs for wheat and many other crops require one or more generations of seed increase before replicated yield trials can be sown. Extensive phenotyping at this stage of the breeding cycle is challenging due to the small plot size and large number of lines under evaluation. Therefore, breeders typically rely on visual selection of small, unreplicated seed increase plots for the promotion of breeding lines to replicated yield trials. With the development of aerial high-throughput phenotyping technologies, breeders now have the ability to rapidly phenotype thousands of breeding lines for traits that may be useful for indirect selection of grain yield. We evaluated early generation material in the irrigated bread wheat (Triticum aestivum L.) breeding program at the International Maize and Wheat Improvement Center to determine if aerial measurements of vegetation indices assessed on small, unreplicated plots were predictive of grain yield. To test this approach, two sets of 1,008 breeding lines were sown both as replicated yield trials and as small, unreplicated plots during two breeding cycles. Vegetation indices collected with an unmanned aerial vehicle in the small plots were observed to be heritable and moderately correlated with grain yield assessed in replicated yield trials. Furthermore, vegetation indices were more predictive of grain yield than univariate genomic selection, while multi-trait genomic selection approaches that combined genomic information with the aerial phenotypes were found to have the highest predictive abilities overall. A related experiment showed that selection approaches for grain yield based on vegetation indices could be more effective than visual selection; however, selection on the vegetation indices alone would have also driven a directional response in phenology due to confounding between those traits. A restricted selection index was proposed for improving grain yield without affecting the distribution of phenology in the breeding population. The results of these experiments provide a promising outlook for the use of aerial high-throughput phenotyping traits to improve selection at the early-generation seed-limited stage of wheat breeding programs.

mapping quality of at least 20 were selected for SNP calling. The SNPs called from the 1 production step were filtered with three criteria: inbreeding coefficient of at least 0.8, 2 Fisher Exact Test (p-value<0.001) to determine biallelic single locus SNPs (Poland et al.,3 2012), and Chi-square test for biallelic segregation with 96 percent expected inbreeding. 4 SNPs that passed at least one filtering criteria were retained, and the subsequent SNP set 5 was filtered to remove those with greater than 50 percent missing data and less than 0.01 6 MAF. Entries with over 50 percent missing data were also removed from the analysis. 7 The final dataset consisted of 13,271 SNP markers for 1,817 entries. Missing markers 8 were imputed with the marker mean. The genomic relationship matrix (G) between 9 entries was calculated according to Endelman and Jannink (2012). The additive pedigree 1 0 matrix (A) was calculated as twice the coefficient of parentage. for each VI at each time-point by fitting the following mixed model: where y ijkl is the trait value; μ is the overall mean; g i is the random genetic effect for the residual effect. For DTHD and DTMT, which were evaluated on the first replicate 1 only, the random effect for replicate was removed from model (2).

2
For the SP and SP BP experiments, BLUPs were calculated for the agronomic traits 3 and VIs as by fitting the model: where y i is the trait value; μ is the overall mean; g i~i id N(0, σ g 2 ) is the random genetic 6 effect for genotype i; and ε i ~iid N(0, σ ε 2 R i ) is the residual variance for genotype i, where 7 R i is the correlation matrix for residual effects parameterized as a one-dimensional 8 autoregressive (AR1) process in the column direction to account for potential spatial 9 correlation of observations. An AR1 process model was not applied in the direction of 1 0 rows as full-sibs were sown adjacent to one another in rows, resulting in confounding 1 1 between genetic and spatial variation in the row direction.

2
To account for the effect of phenology, BLUPs for GY and the VIs in the YT and 1 3 SP were calculated a second time by adding a fixed effect covariate for DTHD to models 1 4 (2) and (3). BLUPs for the VIs in the SP BP were corrected for DTHD including the 1 5 qualitative DTHD scores as a fixed effect covariate in model (3).

6
Significant outliers (p-value < 0.05) were identified using Studentized Residuals 1 7 and were removed from the analysis. In addition, the dataset was subset to remove two 1 8 families that contained disproportionately high numbers of full-sibs. After further 1 9 removing entries based on missing genotypic data, the final dataset included 839 and 920 2 0 entries in the 2016-17 and 2017-18 breeding cycles, respectively.

1
To avoid shrinkage on the same data twice (once during the calculation of iid 2 2 BLUPs and again in the GS models), BLUPs were de-regressed by dividing by their Weights for the error variances of the BLUPs were calculated as follows: where H 2 is the proportion of the total variance explained by the genotypic variance between the YT and SP. Without weights, this information would be otherwise ignored in 9 a two-step genomic prediction procedure. where y is a vector of de-regressed BLUPs calculated from models (2) or (3); X is the 1 5 fixed effects design matrix; b is a vector of fixed effects; Z is the random effects design 1 6 matrix; g is a vector of random effects for genotype estimated assuming g~N(0, Gσ a 2 ) 1 7 where G is the genomic relationship matrix and σ a 2 is the additive genetic variance; and e 1 8 is a vector of residuals where e~iid N(0, Iσ e 2 ), σ e 2 is the residual variance, and I is the 1 9 identity matrix. Narrow-sense heritability was calculated as Trait phenotypic correlations were calculated as the Pearson's correlations 2 1 between the iid BLUPs for each trait derived from models (2)  were used to calculate the genetic correlation between each pair of traits according to: To assess the ability to predict GY of the YT using solely HTP information from 1 0 the SP stage, the following linear regression model was fit using GY and HTP 1 1 information from the YT for the entries in the TRN set: where y i is the de-regressed BLUP of genotype i for GY in the YT calculated from model 1 4 (2), μ is the overall mean, z i is the de-regressed BLUP of genotype i for an VI time-point This model was then applied to predict GY for the breeding lines in the TST set.

8
The de-regressed BLUPs for the same VI time-point observed in the SP calculated from for each VI on each HTP time-point. Univariate Genomic Selection 1 As a basis for comparison, univariate GS models using genomic marker or 2 pedigree information were developed to predict GY for each breeding cycle following the 3 form: where y i are the de-regressed BLUPs for GY in the YT calculated in model (2); μ is the 6 overall mean; g i is the random effect for genotype i with g i~N (0, Gσ g 2 ) or g i~N (0, Aσ g 2 ), 7 where G is the genomic relationship matrix and A is the pedigree matrix; and ε i ~iid N(0,  calculated from model (2); X is the fixed effects design matrix, which is the same for GY and the VI trait. This model was fit assuming

Multi-Trait Genomic Selection
, where G is the genomic relationship matrix, A is the pedigree matrix, Due to the strong associations between GY, VIs, and phenology observed in this 1 4 population, two approaches were evaluated to account for the influence of phenology on 1 5 GY. In the first approach, iid BLUPs for DTHD observed in the YT calculated from 1 6 model (2) were included in prediction models (8), (9), and (5) as a fixed effect covariate.

7
For model validation, the predicted values were correlated to iid BLUPs for GY in the 1 8 YT that had likewise been corrected for DTHD by including DTHD records observed in 1 9 the YT as a fixed effect covariate in model (2). Prediction models (8), (9), and (5) were 2 0 tested both with and without this correction for DTHD.

1
In the second approach, a desired gains selection index taking into account the 2 2 genetic correlation between GY and DTHD was developed to assess the ability of   To validate the prediction models in the context of the restricted selection index, 1 4 models (8), (9), and (5) were fit as described to predict GY for the TST set. The GY 1 5 predictions were then combined with observed records for DTHD in the SP to obtain the 1 6 "predicted" restricted gains index values: for DTHD in the SP for genotype i (Smith, 1936;Hazel, 1943). All phenotypic and genotypic data required to confirm the results presented in this

9
The distributions of BLUPs for GY in the YT were consistent across the two 2 0 breeding cycles, while BLUPs for GY in the SP were more variable within and across 2 1 cycles (Table 1). Both the YT and SP reached the heading and maturity stages around the the SP from 2016-17 but had no effect in the SP from 2017-18. For a given VI time- which had a smaller plot size than the SP and was measured in a different breeding cycle 1 3 than GY in the YT.

4
Correcting for DTHD had a minimal impact on the correlation between NDVI and 1 5 GY when measured both traits were measured in the SP. The effect was more 1 6 pronounced for the YT as well as when NDVI and GY were evaluated in the SP and YT, 1 7 respectively. Strong relationships were also observed between the NDVI time-points and  in the SP would have had no effect on DTHD in the YT.

1
When accounting for DTHD in the NDVI measurements in the SP, the response 1 2 to selection for GY in the YT was reduced to below 1 percent at most levels of selection 1 3 intensity. Accounting for DTHD in the GY records in the SP also reduced the projected 1 4 response to selection for GY in the YT, however the decrease was less pronounced. As  selection. Among the total population of 293, over half were scored as "early" while the 1 remaining were divided relatively evenly among "mid" and "late" scores for DTHD, with 2 less than 2 percent scored as "very late". The distribution for DTHD among the selected 3 lines approximately reflected this, though all "very late" lines were eliminated.

4
In the SP BP , the NDVI time-points had correlations with GY of 0.13 for TP-1 and 5 0.05 for TP-2 (Fig. 3A). Despite these weak relationships, NDVI was effective in 6 identifying SP BP breeding lines with low GY for culling. If the level of selection intensity 7 that the breeders used for visual selection had been applied to the SP BP using the NDVI TP-2, respectively. However, over half of these lines were scored as "late" or "very late" 1 3 for DTHD.

4
Since CIMMYT breeders attempt to avoid selecting directionally on DTHD when 1 5 performing visual selection, a more appropriate comparison would be to identify the top 1 6 breeding lines in terms of NDVI for each qualitative DTHD score. The breeders visually 1 7 selected 22, 8, and 6 lines scored as "early," "mid," and "late" for DTHD, respectively.

8
By taking subsets of the "early," "mid," and "late" lines and identifying the top 22, 8, and 1 9 6 lines, respectively, based on NDVI, the response to selection in GY for the SP BP would 2 0 have been 2.91 percent for TP-1 and 4.23 percent for TP-2 (Fig. 3B). There was still 2 1 minimal overlap between the lines selected by the breeders and the top lines for NDVI 2 2 using this approach. The predictive abilities of the three VIs gave similar results with no individual VI 2 providing a consistent, significant advantage. In addition, the predictive abilities of 3 models using the pedigree relationship matrix were low, with an average of 0.16 for the 4 univariate GS models, likely due to the large numbers of full-sibs and missing 5 information in the pedigree records. Therefore, only the predictive abilities of models 6 using NDVI and the genomic relationship matrix (G) are reported and discussed further. time-points (Fig. 4). This exceeded predictive abilities for univariate GS, though multi-1 1 trait GS models that incorporated both NDVI and genomic information provided the 1 2 highest predictive abilities overall.
1 3 When estimating model predictive abilities, two approaches were taken to account 1 4 for confounding of DTHD with both NDVI and GY. In the first approach, iid BLUPs for 1 5 DTHD in the YT were included as a fixed effect in the prediction model. This correction 1 6 had a negligible impact on univariate GS, but the predictive abilities of the HTP selection 1 7 models were strongly reduced to the extent that they became less predictive than 1 8 univariate GS. Following the correction, the predictive abilities of the multi-trait GS 1 9 models were more similar to those of univariate GS.

0
The second approach was to develop a restricted selection index to identify GS showed the highest predictive abilities for both breeding cycles. The aerial VI traits were found to be more heritable than GY in both the YT and 1 7 SP, and minimal differences were observed with respect to heritability among the knowledge, this is the first report aiming to assess the heritability of aerial VI traits smaller plot sizes than the SP. Together, these results suggest that the VI traits measured 1 7 in at the early-generation seed-limited stages may still be effective for improving GY at result of planting full-sibs adjacent to one another.

0
Harvesting 300 plots of the SP BP for GY enabled a direct comparison between previously evaluated. In this study, GY and HTP records assessed in the YT were 1 3 integrated with genomic information to train multi-trait GS models. These were then 1 4 applied to predict GY for the SP, where collecting HTP traits is possible but GY is not 1 5 typically recorded. Overall, the multi-trait models showed the highest predictive abilities 1 6 as compared to univariate GS and prediction with the HTP traits alone, though the 1 7 increase over the latter was marginal. Therefore, it may be more cost-effective for 1 8 breeding programs to invest in HTP alone as opposed to both HTP and genotyping for 1 9 GS when selecting at the early-generation seed-limited stages. While the cost of breeding programs at the early stages can make applying GS expensive. However, if both 2 2 3 1 genomic markers and HTP traits are available, using both sources of information to build 1 prediction models would be optimal.
2 Multiple findings from this study suggested that using NDVI as a selection variance that may have been contributing to both traits is partitioned into the fixed effect.

0
Selection indices, which take into account genetic correlations between traits, may 2 1 therefore be a more suitable approach when traits are highly related. In this study, visual 2 2 selections made in the SP BP showed that the breeders are simultaneously accounting for increases in GS predictive ability, it is likely that breeders will have access to more that the VI traits could be used in conjunction with qualitative DTHD scores to identify 1 5 breeding lines with higher GY potential for the "early", "mid", and "late" heading 1 6 subgroups. This work was supported through the US Agency for International Development