
==== Front
Genetics
Genetics
genetics
Genetics
0016-6731
1943-2631
Oxford University Press US

38469622
10.1093/genetics/iyae037
iyae037
Investigation
Statistical Genetics and Genomics
AcademicSubjects/SCI01180
AcademicSubjects/SCI01140
Featured
Spatio-temporal modeling of high-throughput multispectral aerial images improves agronomic trait genomic prediction in hybrid maize
https://orcid.org/0000-0002-3600-0657
Morales Nicolas Plant Breeding and Genetics Section, School of Integrative Plant Science, Cornell University, Ithaca, NY 14853, USA

Anche Mahlet T Plant Breeding and Genetics Section, School of Integrative Plant Science, Cornell University, Ithaca, NY 14853, USA

Kaczmar Nicholas S Plant Breeding and Genetics Section, School of Integrative Plant Science, Cornell University, Ithaca, NY 14853, USA

Lepak Nicholas United States Department of Agriculture-Agricultural Research Service, Robert W. Holley Center for Agriculture and Health, Ithaca, NY 14853, USA

https://orcid.org/0000-0003-4537-1182
Ni Pengzun Plant Breeding and Genetics Section, School of Integrative Plant Science, Cornell University, Ithaca, NY 14853, USA
College of Bioscience and Biotechnology, Shenyang Agricultural University, Shenhe District, Shenyang, Liaoning Province, PR China

https://orcid.org/0000-0001-9309-1586
Romay Maria Cinta Institute for Genomic Diversity, Cornell University, Ithaca, NY 14853, USA

https://orcid.org/0000-0002-4351-4023
Santantonio Nicholas Plant Breeding and Genetics Section, School of Integrative Plant Science, Cornell University, Ithaca, NY 14853, USA
School of Plant and Environmental Sciences, Virginia Tech, Blacksburg, VA 24061, USA

https://orcid.org/0000-0002-3100-371X
Buckler Edward S Plant Breeding and Genetics Section, School of Integrative Plant Science, Cornell University, Ithaca, NY 14853, USA
United States Department of Agriculture-Agricultural Research Service, Robert W. Holley Center for Agriculture and Health, Ithaca, NY 14853, USA
Institute for Genomic Diversity, Cornell University, Ithaca, NY 14853, USA

https://orcid.org/0000-0001-6896-8024
Gore Michael A Plant Breeding and Genetics Section, School of Integrative Plant Science, Cornell University, Ithaca, NY 14853, USA

https://orcid.org/0000-0001-8640-1750
Mueller Lukas A Plant Breeding and Genetics Section, School of Integrative Plant Science, Cornell University, Ithaca, NY 14853, USA
Boyce Thompson Institute, Ithaca, NY 14853, USA

https://orcid.org/0000-0001-9522-9585
Robbins Kelly R Plant Breeding and Genetics Section, School of Integrative Plant Science, Cornell University, Ithaca, NY 14853, USA

Sillanpää M Editor
Corresponding author: Plant Breeding and Genetics Section, School of Integrative Plant Science, 102b Beebe Hall, Cornell University, 110 Arboretum Rd, Ithaca, NY 14850, USA. Email: krr73@cornell.edu
Conflicts of interest: The author(s) declare no conflict of interest.

5 2024
12 3 2024
12 3 2024
227 1 iyae03702 12 2023
18 2 2024
24 4 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of The Genetics Society of America.
2024
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Design randomizations and spatial corrections have increased understanding of genotypic, spatial, and residual effects in field experiments, but precisely measuring spatial heterogeneity in the field remains a challenge. To this end, our study evaluated approaches to improve spatial modeling using high-throughput phenotypes (HTP) via unoccupied aerial vehicle (UAV) imagery. The normalized difference vegetation index was measured by a multispectral MicaSense camera and processed using ImageBreed. Contrasting to baseline agronomic trait spatial correction and a baseline multitrait model, a two-stage approach was proposed. Using longitudinal normalized difference vegetation index data, plot level permanent environment effects estimated spatial patterns in the field throughout the growing season. Normalized difference vegetation index permanent environment were separated from additive genetic effects using 2D spline, separable autoregressive models, or random regression models. The Permanent environment were leveraged within agronomic trait genomic best linear unbiased prediction either modeling an empirical covariance for random effects, or by modeling fixed effects as an average of permanent environment across time or split among three growth phases. Modeling approaches were tested using simulation data and Genomes-to-Fields hybrid maize (Zea mays L.) field experiments in 2015, 2017, 2019, and 2020 for grain yield, grain moisture, and ear height. The two-stage approach improved heritability, model fit, and genotypic effect estimation compared to baseline models. Electrical conductance and elevation from a 2019 soil survey significantly improved model fit, while 2D spline permanent environment were most strongly correlated with the soil parameters. Simulation of field effects demonstrated improved specificity for random regression models. In summary, the use of longitudinal normalized difference vegetation index measurements increased experimental accuracy and understanding of field spatio-temporal heterogeneity.

Precisely measuring spatial heterogeneity in field experiments is a persistent challenge in crop research. Here, Morales et al. evaluate approaches to improve spatial modeling using high-throughput phenotypes via unoccupied aerial vehicle imagery. The authors test these approaches with simulation data and data from four Genomes-to-Fields hybrid maize field experiments, finding a two-stage approach improves heritability, model fit, and genotypic effect estimation compared to baseline models.

unoccupied aerial vehicles
spatial correction
genomic prediction
vegetation indices
two-dimensional splines
random regression
autoregressive
permanent environment
high-throughput phenotypes
spatial heterogeneity
soil electrical conductance
elevation
soil curvature
U.S. Department of Agriculture 10.13039/100000199 1024080 100397 1010428 1013637 1013641 Iowa Corn Promotion Board Cornell University 10.13039/100007231
==== Body
pmcIntroduction

Controlling for environmental heterogeneity in agricultural field experiments is critical to obtain accurate estimates of varietal performance and treatment effects (Van Es and Van Es 1993; Brownie et al. 1993; Smith et al. 2005; Xu 2016). In plant breeding where soil composition, elevation, slope, curvature, water content, nutrient availability, and management can vary within field experiments, the genotypic effects driving important agronomic traits can become confounded with a specific plot's permanent environment (PE) effects. In this context, a PE is a plot level nongenetic effect that is persistent across the growing season, giving rise to spatial patterns in the phenotypes of agronomic traits such as yield. Hereafter, PE will be used to describe nongenetic plot level effects estimated using longitudinal data, and spatial effects will refer to nongenetic plot level effects estimated from agronomic data collected only at a single timepoint.

Randomization in experimental designs can help account for confounding genetic, spatial, and environmental variation to a large degree (Piepho et al. 2013; Hoefler et al. 2020), but in early stage trials where replication is limited, it is important to model spatial variation. Statistical approaches, such as the separable autoregressive process and the two-dimensional spline (2DSpl) model, have advanced to capture local dependence effects between experimental plots (Gilmour et al. 1997; Covarrubias-Pazaran 2016). Such spatial effects derive local dependencies from distance-based random covariance structures (e.g. plots that are close to each other are more interdependent than those farther away), but these models often make simplifying assumptions of a consistent rate of decay in interdependency across the entire field. Nonetheless, modeling spatial effects using linear mixed models has improved experimental accuracy in plant breeding (Smith et al. 2005; Robbins et al. 2012; Rodríguez-Álvarez et al. 2018; Bernardeli et al. 2021; Copati et al. 2021).

Spatial heterogeneity can change over the growing season, due to weather and management conditions as well as plant development characteristics. The relationship of time on spatial effects can be explored using models accounting for covariance between timepoints and simple models nested within each timepoint. Repeated measurements in time allow estimation of PE effects from random regression (RR) models, providing a purely temporal representation of spatial heterogeneity. To effectively apply these statistical approaches, measurements should largely span the field; however, it can be difficult and expensive to objectively measure phenotypic traits repeatedly across a field experiment, as seen with quantitative disease resistance traits (Poland and Nelson 2011; Reynolds et al. 2019). Therefore, cost-effective remote sensing approaches are necessary for studying PE in the field across time.

Through the use of unoccupied aerial vehicles (UAVs) and other systems, aerial imaging can reliably and cost-effectively measure high-throughput phenotypes (HTPs) for all experimental plots in the field across the growing season (White et al. 2012; Andrade-Sanchez et al. 2014; Sagan et al. 2019; Sun et al. 2021). A widely studied class of aerial image HTPs are vegetation indices (VI) that include the normalized difference vegetation index (NDVI) (Gitelson et al. 2002; Hunt et al. 2013). VI provide physiologically relevant image features that can track variance such as photosynthetic activity, and have successfully measured chlorophyll content, canopy extent, biomass, and water-use efficiency among other plant attributes (Babar et al. 2006; Bannari et al. 2007; Delegido et al. 2011; Thorp et al. 2018). Promising model-derived HTPs from images exist, such as latent-space and convolutional neural network features; however, our study focused on estimation of a specific plot's NDVI PE from linear mixed models in order to quantify latent spatial heterogeneity and field environmental effects (Taghavi Namin et al. 2018; Gage et al. 2019; Wiesner-Hanks et al. 2019; Feldmann et al. 2021).

While phenotypic data are critical in plant breeding, genomic data are arguably of equal importance. Genomic best linear unbiased prediction (GBLUP) has been extensively applied to predict traits in animals and plants from genome-wide single nucleotide polymorphism (SNP) markers, including in maize and wheat (Meuwissen et al. 2001; Rutkoski et al. 2012; Daetwyler et al. 2013). GBLUP can predict traits from genome-wide SNP marker data by modeling the covariance of additive genetic effects as a genomic relationship matrix (GRM) (VanRaden 2008). Importantly for the models presented in this study, genomic relationships can improve the accuracy of partitioning genetic and nongenetic sources of variation, reducing potential confounding between additive genetic effects and plot level spatial heterogeneity.

VI can improve genomic prediction through multivariate approaches by leveraging genetic correlations between the VI and agronomic traits of interest, as demonstrated for grain yield in wheat and biomass in soybean (Rutkoski et al. 2016; Sakurai et al. 2021). While multitrait models leverage genetic correlations across traits to improve predictions, residual correlations exist between NDVI and grain yield in maize (Anche et al. 2020). Recent studies have successfully proposed two-stage approaches for incorporating HTP, such as detecting spatial effects using the SpATS package in the first stage and then creating P-spline hierarchical growth models in the second stage (Pérez-Valencia et al. 2022). Furthermore, modeling approaches for integrating HTP, genomic information, and environmental information can be generalized into the genotype-to-phenotype model framework (van Eeuwijk et al. 2019); however, modeling NDVI PE for agronomic trait spatial corrections has not previously been widely explored and field tested.

Building on previous work, our study proposed a two-stage approach for improving agronomic trait spatial corrections. To summarize the two-stage approach, the first stage separated NDVI PE from additive genetic effects in the HTP, using either spatial corrections or RR models. The second stage summarized the NDVI PE within GBLUP for the agronomic traits using two distinct implementations, either modeling a plot-to-plot covariance of random effects (L) or fitting regressions on PE estimates from the first stage as fixed effects (FE). The proposed approach studied the following two questions utilizing simulated data and several years of hybrid maize field experiments. Firstly, are NDVI PE consistently able to detect spatial heterogeneity across the growing season affecting end-of-season agronomic traits? Secondly, can NDVI PE be used in the proposed two-stage models to improve spatial corrections for agronomic traits?

Materials and methods

Field experiments

As part of the Genomes-to-Fields (G2F) program, inbred and hybrid maize (Zea mays L.) field evaluations were planted at the Musgrave Research Farm (MRF) in Aurora, NY. Of importance to this study were the hybrid maize experiments planted in 2015 (Genomes to Fields 2020), 2017 (G2F Consortium 2019), 2019 (Genomes to Fields 2021), and 2020 (Genomes to Fields 2022), named 2015_NYH2, 2017_NYH2, 2019_NYH2, and 2020_NYH2, respectively. All experiments were randomized block designs and the tested hybrids varied across experiments (McFarland et al. 2020; AlKhalifah et al. 2018; Lima et al. 2023). The 2015_NYH2, 2017_NYH2, 2019_NYH2, and 2020_NYH2 experiments consisted of 375, 232, 612, and 401 unique hybrids, respectively, encompassing hybrids of diverse genetic backgrounds and expired-proprietary (Ex-PVP) biparental crosses.

The 2015_NYH2 biparental hybrids predominantly featured PHZ51, PHB47, LH82, and LH185 inbred testers and shared a common parent from three inbred biparental cross populations (PHN11 × PHW65, Mo44 × PHW65, and PHW65 × MoG). In 2019_NYH2 many biparental hybrids also shared a common parent from the PHN11 × PHW65 and Mo44 × PHW65 populations; however, PHT69 was the principal tester. The 2017_NYH2 and 2020_NYH2 experiments featured hybrids with shared parents from a MAGIC population (W10004) (Michel et al. 2022). In 2017_NYH2 the same 2015_NYH2 testers were predominantly featured along with 3IIH6 and PHRE1 inbreds; however, in 2020_NYH2, PHP02 was the principal inbred tester. Each experimental plot was seeded in two-row plantings. The 2015_NYH2, 2017_NYH2, 2019_NYH2, and 2020_NYH2 experiments were planted in 100 rows of 10, 10, 16, and 10 ranges, respectively, and plots were evenly distributed across two contiguous blocks in the field. In all field experiments, the following agronomic traits were measured: grain yield (GY) (bu/acre), grain moisture (GM) (%), and ear height (EH) (cm).

Experimental variation across the years arose from the maize hybrids planted, the field sites utilized, the planting and harvest dates, the weather conditions experienced during the growing seasons, and the time points at which imaging events occurred. Supplementary Fig. 1 illustrates the time points at which HTPs were extracted, as well as the growing degree days (GDD) and cumulative precipitation (CP) on those days. GDD were calculated for each imaging event time point from the planting date of the experiment and were found using daily weather data for the ground station at MRF (GHCND:USC00300331) accessed via the NOAA NCEI NCDC database. The evidenced variation in GDD and CP across years highlighted the necessity to collect as many imaging events as possible over the growing season.

Aerial image collection and processing

A MicaSense RedEdge 5-channel multispectral camera mounted onto an UAV captured images in the blue, green, red, near infrared, and red-edge spectra. The UAV flew at an altitude of 25 to 30 m and at a speed of 6 km/h. To complete a flight, the preprogrammed, serpentine flight plans required approximately 35 min to traverse the 3 km path. At least 80% overlap along both axes was ensured in the collected images.

Each flight produced approximately 5,000 images from the MicaSense camera, which were then processed into orthophotomosaics using Pix4dMapper photogrammetry software (Pix4D 2017). To produce reflectance calibrated raster images, the software used the MicaSense radiometric calibration panel images captured immediately prior to each UAV flight, as well as illumination metadata embedded in each capture by the MicaSense camera. Orthophotomosaic, or orthophoto, images were produced with approximately 1 cm per pixel resolution ground sample distance (GSD). The resulting reflectance orthophoto images were then uploaded into ImageBreed software, which enabled plot-polygon templates to be created and assigned to the field experimental design (Morales, Kaczmar, et al. 2020). Supplementary Fig. 2 illustrates a representative near-infrared (NIR) reflectance orthophoto image from 2019_NYH2 taken on August 15, 2019, with the plot-polygons overlaid. NDVI HTP were extracted through ImageBreed, and derived as the mean of NDVI pixel values within the plot boundaries (Gitelson et al. 2002; Hunt et al. 2013; Patrignani and Ochsner 2015; Bhandari et al. 2021). The image, field experiment, phenotypic, and genotypic data within ImageBreed are FAIR and can be queried through openly described APIs (Wilkinson et al. 2016; Selby et al. 2019).

Flights in 2015, 2017, and 2019 were scheduled approximately once per week, while 2020 targeted a frequency of twice per week. Technical problems and poor weather conditions, such as clouds, rain, and high winds, resulted in fewer imaging events being suitable for HTP extraction. To separate the HTP imaging event dates into biological growth stages, GDD ranges were defined to account for the early vegetative phase (P0) at 0 to 1,225 GDD, the active reproductive phase (P1) at 1,226 to 1,800 GDD, and the late reproductive phase (P2) at 1,801 to 2,500 GDD. The three GDD ranges captured periods where NDVI was first steadily increasing, then plateauing, and finally steadily decreasing. Table 1 summarizes the growth stage distribution of imaging event dates for which HTP were successfully extracted.

Table 1. Summary of UAV imaging dates for which HTP were successfully extracted in the 2015, 2017, 2019, and 2020 field experiments.

Field experiment	Planting date	Field location	Growth stage: P0	Growth stage: P1	Growth stage: P2	
2015_NYH2	2015 May 7	MRF, Field P	July 21	Aug 7, Aug 20	Sept 10	
2017_NYH2	2017 May 18	MRF, Field Y	June 12	Aug 2, Aug 17, Sept 1	Sept 6, Sept 12, Sept 24	
2019_NYH2	2019 May 23	MRF, Field N	July 16, July 24, July 29	Aug 5, Aug 15	Sept 10	
2020_NYH2	2020 May 22	MRF, Field D	June 29, July 9, July 15, July 18	July 22, July 28, Aug 1	Aug 20, Aug 26, Sept 9, Sept 18, Oct 2	
Growth stages of P0, P1, and P2 map to early, active, and late reproductive phases broadly defined as 0 to 1,225 GDD, 1,226 to 1,800 GDD, and 1,801 to 2,500 GDD, respectively.

Soil information

In 2019 at the MRF, a ground conductivity meter (EM38-MK2, Geonics, Canada) surveyed Field N where the G2F hybrid maize field experiment was planted. The EM-38 probe used electrical inductance to characterize variation originating from a combination of factors including soil salinity, soil texture, water content, water retention, soil type, and soil nutrients (Heil and Schmidhalter 2017). The goals for the soil information in this study were to (1) better understand driving factors for the estimated NDVI PE from the aerial imagery, and (2) determine whether a soil survey was a practical alternative to aerial imagery for improving agronomic spatial corrections in the second stage. Due to planting rotations and other logistical concerns, the G2F experiments conducted in 2015, 2017, 2019, and 2020 were all in distinct field locations as shown in Table 1; therefore, the soil information in this study could only be applied to the 2019 experiment.

The soil survey was conducted prior to the hybrids being planted by passing the probe over the field in a dual-serpentine pattern, with nine passes in the east–west orientation and 24 passes in the north–south orientation. The georeferenced elevation (Alt) and apparent electrical conductance (EC) data were then interpolated over the entire field using ordinary Kriging in R (Pebesma 2004). The interpolated raster was produced using a spherical variogram model at a resolution of 10−6 WGS84 units covering 120 by 200 cells. Supplementary Fig. 3 illustrates (A) a map of the collected EM38 soil survey, (B) the region to interpolate into, and in (C) and (D) the interpolated soil EC and Alt across the field, respectively. Finally, a mean value for the soil EC and Alt was extracted using ImageBreed for each plot in the 2019_NYH2 field experiment.

In order to approximate soil curvature as elevation gradients, first and second 2D numerical derivatives were computed on the plot level soil EC and altitude measurements, denoted as dEC, d2EC, dAlt, and d2Alt, respectively. 2D numerical derivatives were computed by averaging the differences between a given plot and the three immediately adjacent rings encircling it. Supplementary Fig. 4 illustrated heatmaps of the extracted plot level soil EC and Alt, along with the first and second derivatives.

Genotyping data

G2F maize inbred lines were scored with a genotyping-by-sequencing (GBS) approach that yielded 945,574 SNP markers across the genome for 1,577 samples representing a total of 1,325 unique maize inbred lines (G2F Consortium 2019) (Elshire et al. 2011; McFarland et al. 2020). The resulting genome-wide variant call format (VCF) data were queried for the hybrids in the 2015, 2017, 2019, and 2020 field experiments (Danecek et al. 2011; Morales, Bauchet, et al. 2020; Morales et al. 2022). Due to minor typographical errors (e.g. Mo17 vs MO17), data cleaning was required prior to mapping the sample identifiers in the VCF to the field experiment genotype identifiers and to the pedigree information for the maize hybrids. Genotypes were filtered for SNPs with minor allele frequency <5% or with >40% missing data, and for samples containing >20% missing data. GRMs were computed using the A.mat function in rrBLUP, with missing data imputed as the mean genotype (VanRaden 2008; Endelman 2011). In this study, only additive genetic relationships were modeled (Griffing 1956).

Given that many of the evaluated maize hybrids in the G2F program originated as bi-parental crosses of the genotyped inbred lines and that the inbred pollen and seeds parents are unrelated to each other, hybrid maize genotypes were computed from parental GRMs following:

rij=0.5*(rmi,mj+rpi,pj)

where rij is the genomic relationship between hybrids, rmi,mj is the genomic relationship between the seed parents of the hybrids, and rpi,pj is the genomic relationship between the pollen parents of the hybrids. Therefore, the hybrid genotype is computed by averaging the parents’ genomic relationships for general combining ability (GCA), where the parents’ genomic GCA matrix was calculated as:

GGCA=0.5*(Wt(W)2∑pi(1−pi))

where W is the mean centered allelic dosage matrix (2—homozygous reference allele, 1—heterozygous, 0—homozygous nonreference allele) for the inbred parents and pi is the frequency of the ith SNP marker. If a hybrid's parental inbred lines were not genotyped, then the hybrid was included in the GRM with a diagonal value of one and off-diagonal values of zero.

HTP spatial heterogeneity

The first question of this study was whether a specific plot's nongenetic NDVI PE, which were estimated by spatial and temporal RR effects, consistently detected across the growing season the spatial heterogeneity affecting end-of-season agronomic traits. Therefore, the first stage of the proposed two-stage approach focused on measuring spatial heterogeneity in NDVI across the growing season. All analyses were nested within year given the limited overlap in entries between years.

Variance in the collected HTP observations, VP , was modeled as arising from genetic variance among the maize hybrids, VG , and from environmental sources of variance, VE (Falconer and Mackay 1996). In mixed model matrix notation this was formulated as

(1) y=Xβ+Zaua+Zpup+e

where y is a vector of phenotypic observations and X is an incidence matrix mapping phenotypic values to the fixed effects, β, for replicate (nested within timepoint for analyses including multiple timepoints). The random effects ua and up represented the additive genetic and the specific plot's PE, respectively. The incidence matrices Za and Zp linked the random effects ua and up to the observed phenotypes. The random residual error is represented by e.

The following paragraphs detail how Equation (1) modeled the variation of NDVI PE across time by utilizing either a single HTP time point or many HTP time points simultaneously.

The spatial model fitted either a single HTP time point or multiple HTP time points, a distinction referred to as the single-trait or single-trait-repeated cases. In the single-trait case, spatial models were run independently for each of the collected HTP time points. In contrast, the single-trait-repeated case fitted several collected HTP time points in a single model, such that the vector y represented HTP observations taken at different time points across the growing season:

y=[yt1yt2⋮]

In the single-trait spatial case, the variance of the random additive genetic effect was defined as:

(2) var(ua)=σua2G

where the matrix G represented the GRM between evaluated hybrids (VanRaden 2008). The residual variance was modeled as:

var(e)=σe2I

where I is the identity matrix.

The single-trait-repeated case defined this as

(3) var(ua)=Σua⊗G=[σa12σa1σa2⋯σa2σa1σa22⋯⋮⋮⋱][t,t]⊗G

where the unstructured matrix Σua is of order [t,t] denoting the number of time points involved and ⊗ is the Kronecker product.

var(e)=Σe⊗I=[σe12σe1σe2⋯σe2σe1σe22⋯⋮⋮⋱][t,t]⊗I

where the unstructured matrix Σe is of order [t,t] denoting the number of time points involved, I is an identity matrix with dimensions equal to the number of plots, and ⊗ is the Kronecker product. The advantage of the single-trait-repeated approach was to explicitly model genetic and residual covariances between time points. When model convergence became problematic due to high numbers of time-points in Equation (3), a minimum of three time-points were selected based on NDVI heritability and correlation to yield.

Two methods for modeling environmental spatial variation were investigated, namely 2DSpl and AR1 models. The 2DSpl method used penalized splines fitted using sommer in R (R 3.6.3, Sommer 4.1.3) (Covarrubias-Pazaran 2016). The AR1 method explicitly defined a separable autoregressive covariance structure fitted using ASReml R (R 3.6.3, ASReml R 4.1.0.126) (Gilmour et al. 2015). These two methods use row and column information to model spatial heterogeneity and are further explained in Supplementary File 2. As an alternative to spatial mixed models, random regression (RR) models were explored for separating NDVI PE from additive genetic effects and residual variation specific to a given timepoint (Kirkpatrick et al. 1990; Van der Werf et al. 1998; Schaeffer 2004; Arnold et al. 2019).

RR PE offered a purely longitudinal representation of the plot level spatial heterogeneity in the field. The RR model notably allows estimation of additive genetic effects and PE effects whose covariance follows the GRM and a plot-to-plot correlation matrix E, respectively. The RR model can be expressed in mixed model matrix notation as in Equation (1), however, the incidence matrices Za and Zp contain Legendre polynomial functions of the time, in GDD, at which the measurement was recorded, and ua and up denote the hybrid specific random additive genetic and PE regression coefficients, respectively. The fixed effect was for the design replicate nested with the imaging event date, year, and location, while the random residual e whad heterogeneous variance. The overall variance was written as:

var(y)=Za(Σua⊗G)Z′a+Zp(Σup⊗E)Z′p+Σe⊗I

In this equation, G represents the GRM and E represents a plot-to-plot covariance matrix capturing environmental effects, ideally computed from envirotyping information. Envirotyping aimed to uniquely define the complete environment of an organism by including soil, climate, and developmental parameters; however, this study focused on applying NDVI PE (Xu 2016).

The additive genetic variance for any hybrid at a specific time point (σa,ti2) and the additive genetic covariance for any hybrid between two time points (σa,titj) can be calculated using:

σa,ti2=ztiTΣuazti

and

σa,titj=ztiTΣuaztj

where zti and ztj are vectors of the continuous random regression function evaluated at time points ti and tj, respectively. Similar expressions for the PE variance (σp,ti2) and covariance (σp,titj) can be written as

σp,ti2=ztiTΣupzti

and

σp,titj=ztiTΣupztj

The residual variance structure :

var(e)=Σe⊗I=[σe120⋯0σe22⋯⋮⋮⋱][t,t]⊗I

where the heterogeneous diagonal matrix Σe is of order [t,t] denoting the number of time points involved, I is an identity matrix with dimensions equal to the number of plots, and ⊗ is the Kronecker product. In this study, solutions to the RR model were estimated using the BLUPF90 family of programs (Misztal et al. 2002).

The RR model had computational benefits over the spatial mixed models described previously. Firstly, the number of variance components to estimate was equal to 12n(n−1)+n, where n is the order of the respective random regression functions, regardless of the number of time points t represented in the observations y. Secondly, continuous curves for the random additive genetic and PE effects could be evaluated for any time point because the random regression coefficients fit a covariance function. Thirdly, there was flexibility in the type of regression function that can be fitted, for instance, splines vs Legendre polynomials (Szeg 1975). In this study, third-order Legendre polynomials were considered after looking at the log-likelihood of first- and second-order polynomials as well as the correlation of HTP PE effects and agronomic trait spatial heterogeneity. Furthermore, third-order Legendre polynomials have been found to fit multispectral VI well (Anche et al. 2023). Fourthly, flexible specification of E allowed envirotyping information to be accounted for in the model.

This study explored six structures for the RR PE covariance matrix E. The first approach, named RRID, defined E=I, where I is the identity matrix. The second approach, named RREuc, computed E using the inverse Euclidean distances between experimental plots and standardized between 0 and 1. This is written as

Eij=1(ri−rj)2(ci−cj)2

where Eij is an element of E denoting the relationship between plot i and j, and the variables ri, ci, rj, and cj are the row and column ordinal positions of the plot, respectively. The remaining approaches named RRSoilEC, RRSoilAlt, RR2DSpl, and RRAR1 computed E following:

E=QQ′m

where Q is a centered and standardized matrix of the considered features and m is the number of distinct features. The plot level features in these cases were values of soil EC, dEC, and d2EC, values of soil Alt, dAlt, and d2Alt, values of 2DSplM NDVI random spatial effects, and values of AR1M NDVI random spatial effects, respectively.

The presented models allowed estimation of additive genetic and PE random effects; however, the true genetic and environmental effects were unknown to us. Therefore, to evaluate the robustness of the tested models and their ability to detect field environment features, six different simulation scenarios named Linear, 1D-N, 2D-N, AR1xAR1, and RD were conducted. The simulation methods are described in Supplementary File 1 and were designed to represent purely environmental effects in the field, such as due to soil heterogeneity, altitude, and soil elevation gradients. Each were tested by varying the correlation across time to be 0.75, 0.90, and 1.00 and by setting the simulated variance to be 10%, 20%, and 30% of the total phenotypic variance. Evaluating the accuracy of the simulation process followed:

For a target model, NDVI PE were separated from additive genetic effects present in the real NDVI HTP for a given field experiment.

The computed NDVI PE were subtracted from the NDVI HTP in order to minimize latent spatio-temporal effects in the NDVI.

The target simulation values, meaning one of the six simulation processes, were scaled between 0 and 1, then subsequently scaled to account for either 10%, 20%, or 30% of the observed NDVI phenotypic variation, and finally were added onto the minimized NDVI HTP.

The target model computed PE for the simulation-adjusted NDVI HTP.

Finally, the recovered PE were correlated against the true target simulation values, returning prediction accuracy.

Agronomic spatial correction

The second question in this study was whether the longitudinal NDVI or NDVI PE data could be used to improve spatial corrections for genomic prediction of agronomic traits. Firstly, the following baseline models were defined under the GBLUP framework. The baseline GBLUP model was written as:

(4) y=Xβ+Zaua+e

where y is the agronomic trait of interest, β is a vector of fixed effects for design replication nested by year and location, ua is a vector of the random additive genetic effects of the hybrids, X and Za were incidence matrices linking model terms to y, and e was a vector of the residual errors. The variance of the random additive genetic effect was defined as:

var(ua)=σua2G

where G is the GRM and σua2 is the additive genetic variance. The error variance is defined as:

var(e)=σe2I

where σe2 is the residual variance. Control GBLUP models named G, G + 2DSpl, and G + AR1, are as follows:

(5a) y=XRβR+Zaua+e

(5b) y=XRβR+Zaua+Z2DSplu2DSpl+e

(5c) y=XRβR+Zaua+ZAR1uAR1+e

These models are used as baseline models for agronomic genomic prediction, and defined the vector βR for fixed effects of replication nested by year and location corresponding to incidence matrix XR, while e represented the residual errors. Equation (5a) is the simplest GBLUP case, modeling only random genetic effects ua with corresponding incidence matrix Za. Equations (5b) and (5c) added 2DSpl and AR1 random effects, u2DSpl and uAR1, respectively, with corresponding incidence matrices Z2DSpl and ZAR1. Incorporating spatial corrections accounts for the row and column positions of the plots in the field on which the agronomic trait was measured. Equation (5b) followed the 2DSpl definition from Equation (S1), while Equation (5c) followed the AR1 definition from Equation (S2). Importantly, these three baseline models did not leverage information from the first stage or from the aerial image HTP measurements.

A final baseline defined the multitrait model (M) in Equation (6),

(6) y=XRTβRT+Zaua+e

where y is a vector of both the HTP NDVI time points and the agronomic trait (e.g. GY or EH or GM), βRT is a vector for fixed effects of replicate nested with trait, year, and location, ua is a vector of the random additive genetic effects for all traits, and e is the random residual variance. The random additive genetic (ua) and residual (e) covariances are unstructured across the traits, as in Equation (3). Given the difficulty of fitting large numbers of traits in M, only the two HTP NDVI timepoints with the highest correlation to GY were included in y.

Secondly, to improve on the baseline GBLUP models first-stage, NDVI PE were integrated into the second stage following two implementations. The first implementation modeled the NDVI PE as a plot-to-plot correlation structure for random effects (L), following:

(7) y=XRβR+Zaua+Zpup+e

The variance of the random plot effects up is:

var(up)=σp2L

where:

L=NN′t

The matrix N is the centered and standardized NDVI PE and t is the number of time points considered. The vector βR is for the fixed effects of replicate nested with imaging event date, year, and location, with corresponding incidence matrix XR. This model is similar to Equation (5c) in which the AR1 process explicitly defined a plot-to-plot covariance structure; however, rather than define distance-based assumptions, Equation (7) utilized observed spatial heterogeneity. The second column of Table 2 summarizes the tested L two-stage models in relation to the first-stage.

Table 2. Listed are all nonsoil two-stage models tested in this study.

Stage-1 modela	Stage-2 L model	Stage-2 Havg model	Stage-2 H3 model	Stage-2 Favg model	Stage-2 F3 model	
2DSplU	G + L_2DSplU	G + Havg_2DSplU	G + H3_2DSplU	G + Favg_2DSplU	G + F3_2DSplU	
2DSplM	G + L_2DSplM	G + Havg_2DSplM	G + H3_2DSplM	G + Favg_2DSplM	G + F3_2DSplM	
AR1U	G + L_AR1U	G + Havg_AR1U	G + H3_AR1U	G + Favg_AR1U	G + F3_AR1U	
AR1M	G + L_AR1M	G + Havg_AR1M	G + H3_AR1M	G + Favg_AR1M	G + F3_AR1M	
RRID	G + L_RRID	G + Havg_RRID	G + H3_RRID	G + Favg_RRID	G + F3_RRID	
RREuc	G + L_RREuc	G + Havg_RREuc	G + H3_RREuc	G + Favg_RREuc	G + F3_RREuc	
RR2DSpl	G + L_RR2DSpl	G + Havg_RR2DSpl	G + H3_RR2DSpl	G + Favg_RR2DSpl	G + F3_RR2DSpl	
RRAR1	G + L_RRAR1	G + Havg_RRAR1	G + H3_RRAR1	G + Favg_RRAR1	G + F3_RRAR1	
The first column lists the first-stage models used to separate additive genetic effects from NDVI PE effects. Subsequent columns list models where the second-stage was implemented as a plot-to-plot covariance (L), as a continuous regression (Havg), as three distinct continuous fixed effects (H3), as three distinct binned fixed effects (F3), and as the average of the three binned effects (Favg).

a 2DSplU, two-dimensional spline spatial model fit to single timepoints; 2DSplM, two-dimensional spline spatial model fit to multiple timepoints; AR1U, first-order autoregressive spatial model fit to single timepoints; AR1M, first-order autoregressive spatial model fit to multiple timepoints; RRID, random regression using an identity matrix to model plot-to-plot correlations for PE; RREuc, random regression using Euclidian distance to model plot-to-plot correlations for PE; RR2DSpl, random regression using 2DSpl spatial effect estimates to model plot-to-plot correlations for PE; RRAR1, random regression using first-order autoregressive spatial effect estimates to model plot-to-plot correlations for PE.

Alternatively, fixed FE were defined and followed four variations:

(8a) y=XRβR+β˙HavgHavg+Zaua+e

(8b) y=XRβR+βH˙1H1+βH˙2H2+βH˙3H3+Zaua+e

(8c) y=XRβR+βF˙avgFavg+Zaua+e

(8d) y=XRβR+βF˙1F1+βF˙2F2+βF˙3F3+Zaua+e

In Equation (8a), the vector Havg performed a continuous regression on the average first-stage NDVI PE across all time points, resolving the βH˙avg coefficient. Alternatively, Equation (8b) performed a continuous regression on the vectors H1, H2, and H3 by splitting time points into the previously defined P0, P1, and P2 growth phases, respectively, and then averaging the first-stage NDVI PE within each. This model resolves coefficients βH˙1, βH˙2, and βH˙3 for the P0, P1, and P2 growth phases, respectively. Equations (8c) and (8d) used binned fixed effects derived by reassigning the first-stage NDVI PE to quartile factors (1 = 0–25%, 2 = 26–50%, 3 = 51–75%, 4 = 76–100%), representing poor, marginal, good, and high performing levels. The Favg vector in Equation (8c) was computed by averaging over all time points, resolving the βF˙avg coefficient. The F1, F2, and F3 vectors in Equation (8d) were computed by splitting the time points into the P0, P1, and P2 growth phases, respectively, resolving βF˙1, βF˙2, and βF˙3 coefficients. Columns 3 to 6 of Table 2 summarize the tested FE two-stage model names in relation to the first-stage.

Soil altitude and EC were measured in 2019, enabling Equations (7), (8a) and (8c) to account for soil information rather than first-stage NDVI PE. It was not possible to test Equation (8b) or (8d) with soil information because only a single soil measurement was collected. Representing Equation (7), a model named G + L_RRSoilEC constructed L using soil EC, dEC, and d2EC features, while a model named G + L_RRSoilAlt constructed L using soil Alt, dAlt, and d2Alt features. Using soil measurements in Equation (8a), models named G + Havg_Soil_Alt, G + Havg_Soil_dAlt, G + Havg_Soil_d2Alt, G + Havg_Soil_EC, G + Havg_Soil_dEC, and G + Havg_Soil_d2EC represented Havg derived from the plot level soil Alt, dAlt, d2Alt, EC, dEC, and d2EC, respectively. Whereas for Equation (8c), models named G + Favg_Soil_Alt, G + Favg_Soil_dAlt, G + Favg_Soil_d2Alt, G + Favg_Soil_EC, G + Favg_Soil_dEC, and G + Favg_Soil_d2EC represented Favg derived from the plot level soil Alt, dAlt, d2Alt, EC, dEC, and d2EC, respectively.

To measure the impact of spatial corrections on plot level data points, the plot level heritability is calculated as:

h2=σua2σua2+σe2

where σua2 is the additive genetic variance component and σe2 is the residual error variance component.

Model fit and genotypic effect estimation across replicates (GEER) were used to measure model accuracy. Model fit was defined as the correlation between fitted model predictions y^ and the true phenotypic values y, written as cor(y^,y). Finally, GEER, written as cor(grep1,grep2), involved first partitioning the agronomic trait datasets by replicate, then running the target model on each replicate, and finally correlating the random genetic effect estimates, ua, across the replicates. Given mixed model analyses were performed with varying covariance structures, entry mean heritability was not used as a measure of accuracy.

Results and discussion

Local environmental effects

Before turning attention to the NDVI HTP, spatial corrections for the agronomic traits of GY, GM, and EH were computed. As defined in Equation (5b), Fig. 1 illustrates heatmaps of the 2DSpl spatial effects over the rows and columns of the experimental plots in the 2017_NYH2, 2019_NYH2, and 2020_NYH2 field experiments. Each year the trial was planted in a distinct field location. The heatmaps resolved major, poorly performing regions for GY and EH centered around the row-column positions of (35,2), (19,8), and (25,6) in 2017_NYH2, 2019_NYH2, and 2020_NYH2, respectively. For all three experiments, an inverse spatial pattern was visible for GY and GM, while a similar spatial pattern was visible for GY and EH. The same pattern between traits was evidenced in 2015_NYH2, as illustrated in Supplementary Fig. 7. In 2015_NYH2 a poor performing region for GY and EH was centered near the row-column position of (88,9). In all four years, the proportion of phenotypic variation explained by the 2DSpl spatial effect ranged from ±25 bu/acre of GY, ±2% of GM, and ±15 cm of EH. These results illustrated the importance of spatial heterogeneity on the agronomic traits.

Fig. 1. 2DSpl spatial random effects detected in the agronomic traits of GY, GM, and EH in the 2017, 2019, and 2020 field experiments. The 2015 experiment is shown in Supplementary Fig. 7.

Similar patterns of spatial random effects were found using the 2DSpl and AR1 models defined in Equations (5b) and (5c), respectively. Table 3 lists the correlations between the 2DSpl and AR1 spatial effects for GY, GM, and EH individually in the 2017_NYH2, 2019_NYH2, and 2020_NYH2 experiments. The strongest average correlation was 0.88 for EH, followed by 0.87 for GY, and 0.49 for GM. The 2015_NYH2 experiment was not included in Table 3 because the AR1 model did not converge.

Table 3. Correlations between 2DSpl and first-order autoregressive (AR1) spatial effects for GY, GM, and EH individually in the 2017, 2019, and 2020 field experiments.

Experiment	Grain yield (GY)	Grain moisture (GM)	Ear height (EH)	
2017_NYH2	0.90	0.49	0.74	
2019_NYH2	0.91	0.59	0.93	
2020_NYH2	0.79	0.40	0.97	

To understand correspondence between PE affecting NDVI and the end-of-season agronomic traits, first-stage NDVI PE were compared to GY, GM, and EH spatial effects. Figure 2 shows correlations between the 2DSpl GY spatial effects and the 2DSplU NDVI PE across 12 time points in the 2020_NYH2 field experiment. Illustrated in Fig. 2, correlations at 38 and 48 days after planting (DAP) were 0.6 and 0.7, respectively, and correlations fluctuated between 0.5 and 0.7 throughout the season. Corresponding heatmaps in Fig. 2 illustrate spatial distributions over the rows and columns of the experimental plots, and consistently reveal a large region near the center of the field negatively impacting both GY and NDVI. For reference, phenotypic correlations between NDVI and GY in 2020_NYH2 ranged from a low of 0.04 at 133 DAP to a high of 0.39 at 54 DAP, which were weaker than the correlations observed in Fig. 2. Similarly, the RR model PE effects correlated with the GY spatial effects more strongly than the NDVI and GY themselves. Supplementary Fig. 8 illustrates RRID PE effects correlated against GY and the GY 2DSpl and AR1 spatial effects. The RRID PE tended to be strongly correlated through time and identified similar spatial patterns as the spatial effect models.

Fig. 2. 2020_NYH2 single timepoint 2DSplU NDVI PE estimates observed over 12 time points correlated with GY and GY 2DSpl spatial effects. Corresponding heatmaps showed values over the rows and columns of all experimental plots in the field and revealed similar spatial patterns. As indicated by (1), correlations of 0.5 to 0.7 between GY 2DSpl and the 2DSplU NDVI PE were found throughout the growing season, even early on at 38 DAP.

In 2019_NYH2, 2015_NYH2, and 2017_NYH2 similar correlations between 2DSplU and GY spatial effects were observed, as illustrated in Fig. 3, Supplementary Figs. 9, and 10, respectively. For reference, phenotypic correlations between NDVI and GY in 2015_NYH2, 2017_NYH2 and 2019_NYH2 ranged from a low of 0.13 at 126 DAP, 0.17 at 25 DAP, and 0.33 at 110 DAP, respectively, to a high of 0.42 at 92 DAP, 0.49 at 91 DAP, and 0.67 at 84 DAP, respectively. Therefore, in all experiments the NDVI PE were more correlated to the GY 2DSpl effects across the growing season than the NDVI were correlated to GY. Furthermore, in all years, the spatial patterns affecting GY and EH were detectable by NDVI PE to a large degree (>0.5 correlation), evidenced early on in the growing season at 76 DAP or less.

Fig. 3. 2019_NYH2 single timepoint 2DSplU NDVI PE estimates observed over six time points correlated with GY and GY 2DSpl spatial effects. Corresponding heatmaps showed values over the rows and columns of all experimental plots in the field and revealed similar spatial patterns. Illustrated are (1) correlations between GY 2DSpl and the 2DSplU NDVI PE, (2) correlations between soil EC and NDVI PE. and (3) correlations between soil d2Alt and NDVI PE. The AR1 analog is found in Supplementary Fig. 11.

As a potential alternative to aerial imaging and to better understand the observed spatial effects, soil EC, Alt, and the first and second 2D numerical derivatives (dEC, d2EC, dAlt, d2Alt) were compared with the NDVI PE in the 2019_NYH2 field experiment. Figure 3 includes correlations and heatmaps of the soil measurements with the spatial effects and demonstrated: (1) the 2DSpl spatial effects of GY in 2019_NYH2 correlated up to 0.7 with 2DSplU NDVI PE, (2) soil EC correlated up to 0.5 with 2DSplU NDVI PE, and (3) soil d2Alt correlated up to 0.3 with 2DSplU NDVI PE. In Fig. 3, the soil elevation gradients, represented by heatmaps of Alt, EC, and their derivatives, highlighted the contours of the observed 2DSpl spatial effects. Supplementary Fig. 11 presents an analog to Fig. 3 showing AR1U NDVI PE, GY, and GY AR1 spatial effects in 2019_NYH2, and showed (1) overall weaker correlations between the GY AR1 spatial effects and the AR1U NDVI PE, with a high of 0.5, (2) similar correlations to soil EC with a high of 0.5, and (3) overall weaker correlations to soil d2Alt with a high of 0.2. Therefore, the AR1U NDVI PE were less correlated with the soil parameters than the 2DSplU NDVI PE, and the AR1U NDVI PE were less correlated with GY AR1 effects than the 2DSplU NDVI PE were correlated with the GY 2DSpl effects.

Table 4 summarizes correlations between the soil information in 2019_NYH2 and the model NDVI PE, illustrating how the average first-stage NDVI PE across all time points correlated to soil EC, dEC, d2EC, Alt, dAlt, and d2Alt. Table 4 indicates 2DSplU had the strongest correlation to EC of 0.46 and also correlated relatively strongly with d2Alt. The d2Alt tended to correlate more strongly than Alt or dAlt with the NDVI PE, indicating the importance of elevation gradients in the field. The correlations between 2DSpl NDVI PE and soil parameters indicated that soil information was capturing similar spatial information as the NDVI aerial imaging.

Table 4. Correlations between average NDVI PE in 2019_NYH2 computed using the listed first-stage models and the soil EC, dEC, d2EC, Alt, dAlt, and d2Alt measurements.

Modela	Soil EC	Soil dEC	Soil d2EC	Soil Alt	Soil dAlt	Soil d2Alt	
2DSplU	0.46	0.21	0.04	−0.06	0.11	0.21	
2DSplM	0.21	0.12	0.06	0.08	0.34	0.29	
AR1U	0.39	0.18	0.04	−0.06	0.06	0.16	
AR1M	0.01	0.04	0.06	0.10	0.25	0.21	
RRID	0.10	0.06	0.04	0.04	0.14	0.17	
RREuc	−0.06	0.03	0.09	0.07	0.14	0.15	
RRAR1	−0.18	−0.06	0.03	0.11	0.23	0.21	
RR2DSpl	−0.19	−0.07	0.02	0.13	0.24	0.23	
RRSoilEC	−0.22	−0.08	0.03	0.12	0.25	0.22	
RRSoilAlt	−0.21	−0.09	0.04	0.12	0.25	0.23	
a 2DSplU, two-dimensional spline spatial model fit to single timepoints; 2DSplM, two-dimensional spline spatial model fit to multiple timepoints; AR1U, first-order autoregressive spatial model fit to single timepoints; AR1M, first-order autoregressive spatial model fit to multiple timepoints; RRID, random regression using an identity matrix to model plot-to-plot correlations for PE; RREuc, random regression using Euclidian distance to model plot-to-plot correlations for PE; RR2DSpl, random regression using 2DSpl spatial effect estimates to model plot-to-plot correlations for PE; RRAR1, random regression using first-order autoregressive spatial effect estimates to model plot-to-plot correlations for PE; RRSoilEC, random regression using soil electrical conductance measurements to model plot-to-plot correlations for PE; RRSoilAlt, random regression using altitude measurements to model plot-to-plot correlations for PE.

Summarizing results between the agronomic trait 2DSpl effects and the tested first-stage model PE, Fig. 4 presents correlations for all agronomic traits and years. Figure 4 indicates the 2DSplU, 2DSplM, and AR1U models produced NDVI PE most correlated to the 2DSpl effects of GY and EH in all years and in nearly all time points, while the RR2DSpl, RRAR1, and RRSoilEC models were most correlated to the 2DSpl effects of GM in 2015, 2017, and 2019. The 2DSplU and AR1U NDVI PE correlated with GY and EH 2DSpl spatial effects greater than 0.5 in all years at 80 to 90 DAP. Significant similarities were seen between models run on traits within a given year, for instance both GY and EH in 2017 showed a large peak at 110 DAP and in 2020 both showed a continuous gradual decline. There was an inverted behavior between the GY and GM spatial effects in all years, describable as: in 2015 a high for GY and a low for GM at 105 DAP, in 2017 a high for GY and a low for GM at 110 DAP, in 2019 a high for GY and a low for GM around 90 to 100 DAP, and in 2020 a gradual decline for GY and a gradual incline for GM. In contrast, there was a similar behavior between the GY and EH spatial effects in all years. Supplementary Fig. 12 illustrates the AR1 analog of Fig. 4 with agronomic trait AR1 spatial effects instead of 2DSpl effects and demonstrated weaker correlations to the NDVI PE in all traits and all years. Highly similar patterns in the correlation curves were observed; however, there was a tendency for the GM AR1 spatial effects to correlate with RRID, RR2DSpl, and RRAR1 NDVI PE more strongly than the 2DSpl or AR1 NDVI PE.

Fig. 4. Correlations between 2DSpl spatial effects of agronomic traits and PE from the tested first-stage models. Spatial effects for GY, GM, and EH in the 2015_NYH2, 2017_NYH2, 2019_NYH2, and 2020_NYH2 field experiments were compared with model PE across the growing season. RR models including soil information rather than NDVI in the first-stage, named RRSoilAlt and RRSoilEC, were also included for 2019_NYH2. The AR1 analog is found in Supplementary Fig. 12.

Simulation tested the efficacy of the first stage in detecting known environmental field effects. In each of the six simulation processes (linear, 1D-N, 2D-N, AR1xAR1, random, and RD) 10 iterations were performed, each time generating a new simulation. The simulated environmental variance was tested at 0.1, 0.2, and 0.3 times the proportion of phenotypic variation, and the correlation between time points was tested at 0.75, 0.90, and 1. Supplementary Figs. 22–24 illustrate the results for all simulation scenarios using 2017_NYH2, 2019_NYH2, and 2020_NYH2 NDVI phenotypes, respectively, demonstrating the impacts of varying the simulation environmental variance as well as the correlation of simulated environmental effects across the growing season. Increasing the variance tended to slightly increase the prediction accuracy, while decreasing the correlation between time points tended to decrease prediction accuracy. Figure 5 aggregates the prediction accuracies for the linear, 1D-N, 2D-N, AR1xAR1, and RD simulation scenarios and illustrated model groupings determined by a Tukey Honest Significant Difference (HSD) test.

Fig. 5. Aggregated prediction accuracy of the tested first-stage models for the linear, 1D-N, 2D-N, AR1xAR1, and RD simulation scenarios using the 2017_NYH2, 2019_NYH2, and 2020_NYH2 NDVI data. Prediction accuracy is the correlation of the simulated environmental effect and the model's recovered environmental effect. Models are grouped together after performing a Tukey's HSD test.

The RREuc model showed relatively poor performance. This is possibly due to a mismatch in the geometry of the experiment because in reality the plots in the field were rectangular (e.g. 10 ft by 3 ft) and not perfectly square. The RREuc model was also most sensitive to the tested years and to changes in simulated variance and correlation. The AR1U model performed best in the AR1xAR1 scenario, while the 2DSplU model performed well in the Linear scenario. Specifying the PE covariance matrix allowed the RRAR1 and RR2DSpl models to perform consistently well in the 1D-N and 2D-N scenarios; however, by making no assumptions on the spatial structure the RRID model performs on average less than 10% worse and with comparable consistency.

Estimation of spatial effects

First-stage NDVI PE were incorporated into the second stage of the proposed spatial correction approach using two distinct implementations, either modeling L or FE. Supplementary Fig. 13 illustrates genomic heritability, model fit, and GEER in the four years for GY, GM, and EH for all two-stage models when modeling L random effects. Baseline and spatially corrected GBLUP models, named G, G + 2DSpl, and G + AR1 representing Equations (5a), (5b), and (5c), respectively, are shown. Statistical significance of the spatial corrections and two-stage models were compared to the baseline G model using a paired t-test. The best models were determined by ranking the t-test P-value of GEER for GY. The best five two-stage L models, representing Equation (7), were G + L_AR1U, G + L_AR1M, G + L_2DSplU, G + L_2DSplM, and G + L_RRID.

Alternatively, modeling NDVI PE as FE followed four distinct definitions in Equations (8a–8d). Supplementary Fig. 17 illustrates genomic heritability, model fit, and GEER in the four years for GY, GM, and EH for all two-stage models when modeling FE. The baseline and spatially corrected GBLUP models, named G, G + 2DSpl, and G + AR1, respectively, are shown. As before, the best two-stage FE models were determined by ranking the t-test P-value of GEER for GY. The top eight models were named G + Havg_AR1U, G + Havg_2DSplU, G + Havg_2DSplM, G + Favg_2DSplU, G +H3_2DSplM, G + H3_RRID, G + F3_AR1U, and G + F3_2DSplU.

To measure performance of the proposed two-stage approach over the baseline G model, a difference (G Diff) was computed within each of the four years for heritability, model fit, and GEER. In the following analyses, the baseline G and G + 2DSpl models include results from 2015, 2017, 2019, and 2020, while the baseline G + AR1 model excludes 2015 due to AR1 convergence issues. Figure 6 illustrates G Diff for the baseline G + 2DSpl and G + AR1 spatial correction models and for the best six models defining L (G + L) and FE (G + H). Supplementary Figs. 14 and 18 illustrate all two-stage models when defining L and FE, respectively. The spatially corrected baseline models, G + 2DSpl and G + AR1, demonstrate improvements over G. The G + 2DSpl model provided significant improvements in heritability and model fit for all traits, and a significant improvement in GEER for GY; while, G + AR1 provided significant improvements in model fit for all traits, and a significant improvement in heritability and GEER for GY.

Fig. 6. Differences compared to G (G Diff) for genomic heritability, model fit, and GEER for GY, GM, and EH in the four years. The models G, G + 2DSpl, and G + AR1 were baseline GBLUP and spatially corrected GBLUP models using 2DSpl and first-order autoregressive spatial models, respectively. M is a baseline multitrait model. Illustrated are the best three two-stage models using L (G + L) and the best three two-stage models using FE (G + H), determined by ranking t-test P-value of GEER for GY. The G + L and G + H models have L and FE, respectively, defined using NDVI PE of corresponding names.

Figure 6 demonstrates further improvements for the two-stage models over the baseline G, G + 2DSpl, and G + AR1 models. Two-stage models incorporating NDVI PE improved heritability and GEER for GY, GM, and EH more than the baseline spatially corrected models. The best two-stage models translated increased heritability to an increase in GEER. Improvements to GY GEER over baseline G, G + 2DSpl, and G + AR1 models were summarized in Table 5 for the six best two-stage models when incorporating NDVI PE. Table 5 columns of “ΔG”, “ΔG + 2Dspl”, and “ΔG + AR1” indicate the mean and standard deviations of model differences in GEER for GY for the four years compared to G (G Diff), G + 2DSpl (G + 2DSpl Diff), and G + AR1 (G + AR1 Diff), respectively. Supplementary Figs. 15 and 19 illustrate the G + 2DSpl Diff for all models when defining L and FE, respectively. Supplementary Figs. 16 and 20 illustrate the G + AR1 Diff for all models when defining L and FE, respectively. While GEER for GY and EH was improved, none of the FE models significantly improved GEER for GM, a result potentially attributable to the lower correlations between NDVI PE and GM spatial effects seen in Fig. 4 and Supplementary Fig. 12, the difficulty in detecting GM spatial effects seen in Table 3, and the small GM spatial variation seen in Fig. 1.

Table 5. The best two-stage models vs the baseline GBLUP models (G, G + 2DSpl, G + AR1) when comparing the correlation of GEER for GY.

Modela	ΔG	ΔG + 2DSpl	ΔG + AR1	
G + H3_RRID	0.188 ± 0.094 (*)	0.095 ± 0.045 (*)	0.082 ± 0.053 (+)	
G + H3_2DSplM	0.123 ± 0.055 (*)	0.029 ± 0.074	0.025 ± 0.09	
G + Havg_2DSplU	0.12 ± 0.056 (*)	0.026 ± 0.048	0.012 ± 0.057	
G + L_2DSplM	0.065 ± 0.049 (*)	−0.028 ± 0.105	−0.056 ± 0.123	
G + L_AR1U	0.132 ± 0.077 (*)	0.038 ± 0.035 (+)	0.041 ± 0.031 (+)	
G + L_RRID	0.133 ± 0.041 (*)	0.039 ± 0.043 (+)	0.036 ± 0.049	
The symbols (*) and (+) denote t-test P-values less than 0.05 and 0.1, respectively.

a H3, second-stage spatial corrections modeled as fixed regressions on first-stage spatial effect estimates from three distinct growth phases; Havg, second-stage spatial corrections modeled as a fixed regression on the average first-stage spatial effect estimates; L, second-stage spatial corrections modeled using a plot-to-plot correlation matrix calculated using first-stage spatial effect estimates; RRID, first-stage random regression model using an identity matrix to model plot-to-plot correlations; 2DSplM, first-stage, multitimepoint model using 2DSpl to estimate spatial effects; 2DSplU, first-stage, single time point model using 2DSpl to estimate spatial effects; AR1U, first-stage, single time point model using first-order autoregressive correlations to estimate spatial effects.

Drawing from the simulation results in Fig. 5, the 2DSplM and AR1M models may have had less overall accuracy due to restricted estimation of spatial covariance components between time points, and thereby negatively impacted the simulation when weaker correlations (< 0.9) across time points were used. In real data, as seen in Figs. 2, 3, Supplementary Figs. 8, 9, and 10, NDVI PE tended to be strongly correlated (>0.8) between time points; this may have explained the improvements seen in Table 5 when the 2DSplM NDVI PE were incorporated into second-stage genomic prediction. The RRAR1 and RR2DSpl models performed well in first-stage simulation; however, the second-stage genomic prediction was not particularly improved by these models. The RRID, AR1U, and 2DSplU models performed well in first-stage simulation and significantly improved the second-stage model accuracy, indicating these models provided robust detection of spatial heterogeneity.

The second stage in the proposed approach could use soil data as an alternative to NDVI PE from the first-stage. Figure 7 illustrates differences in 2019_NYH2 against the baseline G (G Diff) for genomic heritability, model fit, and GEER for GY, GM, and EH. The best eight models when modeling L (G + L) or FE (G + H) using soil data are illustrated in Fig. 6; however, Supplementary Fig. 21 illustrates all of the models using soil data. Again, the baseline spatially corrected models, G + 2DSpl and G + AR1, are shown. Included were the L models named G + L_RRSoilEC and G + L_RRSoilAlt, and the FE models named G + Favg_Soil_Alt, G + Favg_Soil_dAlt, G + Favg_Soil_d2Alt, G + Havg_Soil_EC, G + Havg_Soil_dEC, and G + Havg_Soil_d2EC. The soil information increased GM and EH heritability, and model fit for all traits, more than the baseline G + 2DSpl and G + AR1 models; however, for all traits GEER performed lower than the baseline G + 2DSpl and G + AR1 models, particularly for GM and EH.

Fig. 7. Differences compared to G (G Diff) for genomic heritability, model fit, and GEER in 2019_NYH2 for GY, GM, and EH, with soil data implemented as a plot-to-plot correlation matrix (g + l) or as fixed effects (g + h). The best eight models using soil altitude (Alt), soil electrical conductance (EC), and the first and second derivatives (dAlt, d2Alt, dEC, d2EC) were illustrated.

This approach to incorporate soil data did not improve model genotypic effect accuracy for GY; however, this result was from a single field experiment in a single year. Similar to how the NDVI HTP itself did not correlate highly with GY while the spatial effects of NDVI correlated strongly with the spatial effects of GY, the spatial effects of the soil data may prove more beneficial for improving prediction accuracy than the soil data itself. The soil data had much weaker correlations than the NDVI to agronomic traits, with a high of 0.07, 0.03, and 0.20 for GY, GM, and EH, respectively, compared to NDVI with a high of 0.67, 0.60, and 0.46 for GY, GM, and EH, respectively. Further indicating persistent spatial effects may be limiting the effectiveness of soil EC data in this study, Table 4 illustrates that including the soil data into RR models resulted in NDVI PE relatively well correlated with the soil elevation gradients, but negatively correlated with the soil EC data itself. Furthermore, soil data may need to be incorporated with weather information in order to effectively estimate the benefit or detriment of the local environmental effect. For instance, low elevation can be either beneficial or detrimental depending on rainfall.

Conclusion

In all years and for all agronomic traits, correlations between the agronomic trait spatial effects and NDVI PE were higher than correlations between the agronomic traits and NDVI themselves, indicating the spatial patterns in NDVI do provide information on spatial patterns observed for key agronomic traits. Furthermore, the NDVI PE from 2DSpl, AR1, and RR models consistently identified the same poorly performing regions in the field over the growing season, and identified substantially the same regions as the baseline GY and EH spatial effects. Baseline GM spatial effects showed an inverted behavior with NDVI PE and were less localized than for GY and EH. The soil EC correlated most with NDVI PE from the 2DSplU and AR1U models across time. Therefore, spatial heterogeneity quantified by NDVI PE corresponded strongly with agronomic trait spatial effects and soil EC.

Incorporating first-stage NDVI PE into the second-stage spatial corrections for GY, EH, and GM either as a covariance of random effects (L) or as FE, significantly improved heritability, model fit, and GEER. In simulation, the RRAR1 and RR2DSpl models performed strongly; however, only the RRID, AR1U, and 2DSplU models performed well in simulation and also improved the two-stage spatial correction for agronomic traits. The RRID model performed consistently above average in simulation and improved spatial correction performance of GY and EH experimentally. Notably, the RRID model made no spatial assumptions, which may facilitate deployment compared to models which rely on spatial annotations. The equilibrium between model generalizability and model over-specification when detecting NDVI PE was balanced most by the RRID, AR1U, and 2DSplU models.

Aerial image HTP provided greater understanding of spatial heterogeneity in the field, and when coupled into the proposed two-stage spatial correction approach, enabled a more effective spatial correction than any of the baseline models (G + 2DSpl, G + AR1, and M). Furthermore, the observed spatial heterogeneity could be partially explained using soil EC and elevation. Further research is needed to identify more informative image features and develop novel statistical approaches for integrating HTP across the growing season with end-of-season agronomic trait prediction. To these ends, larger datasets are required to evaluate the proposed approaches, and the continued aggregation of FAIR data is crucial.

Supplementary Material

iyae037_Supplementary_Data

iyae037_Peer_Review_History

Acknowledgments

The authors would like to thank the G2F consortium for providing the NYH2 field experiments from 2015 to 2020 used in this study. This consortium involves more than 30 researchers representing more than 20 research institutions. Details about the initiative and publicly available resources can be found at www.Genomes2Fields.org. Thanks to Chris Hernandez, Peter Selby, Sam Bouabane, Ranjita Thapa, Simon Reinhard, Lynn Johnson, Seth Murray, Jacob Washburn, Filipe I. Matias, Annarita Marrano, and Felipe Sabadin for their help and suggestions on the image processing pipeline and the research more broadly.

Data availability

This study used phenotypic data of hybrid maize (Z. mays L.) field experiments part of the G2F program planted in 2015 (https://doi.org/10.25739/erxg-yn49), 2017 (https://doi.org/10.25739/w560-2114), 2019 (https://doi.org/10.25739/t651-yy97), and 2020 (https://doi.org/10.25739/hzzs-a865), named 2015_NYH2, 2017_NYH2, 2019_NYH2, and 2020_NYH2, respectively. The genotypic SNP marker data was also from the G2F program (https://doi.org/10.25739/frmv-wj25). The collected image data from 2015, 2017, 2019, and 2020 are available in the Supplemental section of this manuscript.

Supplemental material available at GENETICS online.

Funding

This work was supported by the U.S. Department of Agriculture National Institute of Food and Agriculture, Hatch projects 1024080 (K.R.R.), 100397 (M.A.G), 1010428 (M.A.G.), 1013637 (M.A.G.). 1013641 (M.A.G.), Iowa Corn Promotion Board, and Cornell University startup funds (K.R.R. and M.A.G.).
==== Refs
Literature cited

AlKhalifah  N, Campbell  DA, Falcon  CM, Gardiner  JM, Miller  ND, Romay  MC, Walls  R, Walton  R, Yeh  CT, Bohn  M, et al  2018. Maize genomes to fields: 2014 and 2015 field season genotype, phenotype, environment, and inbred ear image datasets. BMC Res Notes. 11 (1 ):452. doi:10.1186/s13104-018-3508-1.29986751
Anche  MT, Kaczmar  NS, Morales  N, Clohessy  JW, Ilut  DC, Gore  MA, Robbins  KR. 2020. Temporal covariance structure of multi-spectral phenotypes and their predictive ability for end-of-season traits in maize. Theor Appl Genet. 133 (10 ):2853–2868. doi:10.1007/s00122-020-03637-6.32613265
Anche  MT, Morales  N, Kaczmar  NS, Santantonio  N, Gore  MA, Robbins  KR. 2023. Scalable growth models for time-series multispectral images. Plant Phenome J. 6 (1 ):e20064. doi:10.1002/ppj2.20064.
Andrade-Sanchez  P, Gore  MA, Heun  JT, Thorp  KR, Elizabete Carmo-Silva  A, French  AN, Salvucci  ME, White  JW. 2014. Development and evaluation of a field-based high-throughput phenotyping platform. Funct Plant Biol. 41 (1 ):68–79. doi:10.1071/FP13126.
Arnold  PA, Kruuk  LEB, Nicotra  AB. 2019. How to analyse plant phenotypic plasticity in response to a changing climate. New Phytol. 222 (3 ):1235–1241. doi:10.1111/nph.15656.30632169
Babar  MA, Reynolds  MP, van Ginkel  M, Klatt  AR, Raun  WR, Stone  ML. 2006. Spectral reflectance to estimate genetic variation for in-season biomass, leaf chlorophyll, and canopy temperature in wheat. Crop Sci. 46 (3 ):1046–1057. doi:10.2135/cropsci2005.0211.
Bannari  A, Shahid Khurshid  K, Staenz  K, Schwarz  JW. 2007. A comparison of hyperspectral chlorophyll indices for wheat crop chlorophyll content estimation using laboratory reflectance measurements. IEEE Trans Geosci Remote Sens. 45 (10 ):3063–3074. doi:10.1109/TGRS.2007.897429.
Bernardeli  A, Rocha  JR, Borém  A, Lorenzoni  R, Aguiar  R, Silva  JN, Bueno  RD, Alves  RS, Jarquin  D, Ribeiro  C, et al  2021. Modeling spatial trends and enhancing genetic selection: an approach to soybean seed composition breeding. Crop Sci. 61 (2 ):976–988. doi:10.1002/csc2.20364.
Bhandari  M, Baker  S, Rudd  JC, Ibrahim  AMH, Chang  A, Xue  Q, Jung  J, Landivar  J, Auvermann  B. 2021. Assessing the effect of drought on winter wheat growth using unmanned aerial system (UAS)-based phenotyping. Remote Sens. 13 (6 ):1144. doi:10.3390/rs13061144.
Brownie  C, Bowman  DT, Burton  JW. 1993. Estimating spatial variation in analysis of data from yield trials: a comparison of methods. Agron J. 85 (6 ):1244–1253. doi:10.2134/agronj1993.00021962008500060028x.
Copati  MG, Dariva  FD, Dias  FD, Rocha  JR, Pessoa  HP, de Almeida  GQ, Carneiro  PC, Nick  C. 2021. Spatial modeling increases accuracy of selection for phytophthora infestans -resistant tomato genotypes. Crop Sci. 61 (6 ):3919–3930. doi:10.1002/csc2.20584.
Covarrubias-Pazaran  G . 2016. Genome-assisted prediction of quantitative traits using the R package sommer. PLoS One. 11 (6 ):e0156744. doi:10.1371/journal.pone.0156744.27271781
Daetwyler  HD, Calus  MPL, Pong-Wong  R, de Los Campos  G, Hickey  JM. 2013. Genomic prediction in animals and plants: simulation of data, validation, reporting, and benchmarking. Genetics. 193 (2 ):347–365. doi:10.1534/genetics.112.147983.23222650
Danecek  P, Auton  A, Abecasis  G, Albers  CA, Banks  E, DePristo  MA, Handsaker  RE, Lunter  G, Marth  GT, Sherry  ST, et al  2011. The variant call format and VCFtools. Bioinformatics. 27 (15 ):2156–2158. doi:10.1093/bioinformatics/btr330.21653522
Delegido  J, Verrelst  J, Alonso  L, Moreno  J. 2011. Evaluation of sentinel-2 red-edge bands for empirical estimation of green LAI and chlorophyll content. Sensors. 11 (7 ):7063–7081. doi:10.3390/s110707063.22164004
Elshire  RJ, Glaubitz  JC, Sun  Q, Poland  JA, Kawamoto  K, Buckler  ES, Mitchell  SE. 2011. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS One. 6 (5 ):e19379. doi:10.1371/journal.pone.0019379.21573248
Endelman  JB . 2011. Ridge regression and other kernels for genomic selection with R package rrBLUP. Plant Genome. 4 (3 ):250–255. doi:10.3835/plantgenome2011.08.0024.
Falconer  DS, Mackay  TFC. 1996. Introduction to Quantitative Genetics. 4th ed. Harlow, UK: Logmans Green.
Feldmann  MJ, Gage  JL, Turner-Hissong  SD, Ubbens  JR. 2021. Images carried before the fire: the power, promise, and responsibility of latent phenotyping in plants. Plant Phenome J. 4 (1 ):e20023. doi:10.1002/ppj2.20023.
G2F Consortium. 2019. GenomesToFields 2014-2017 Datasets, CyVerse Data Commons CyVerse Data Commons. doi:10.25739/frmv-wj25.
Gage  JL, Richards  E, Lepak  N. 2019. In-field whole-plant maize architecture characterized by subcanopy rovers and latent space phenotyping. Plant Phenome. 2 (1 ):1–11. doi:10.2135/tppj2019.07.0011.
Genomes to Fields. 2020. Genomes_to_Field_Planting_Season_2015_V2. CyVerse Data Commons. doi:10.25739/erxg-yn49.
Genomes to Fields. 2021. Genomes to Fields 2019 dataset. CyVerse Data Commons. doi:10.25739/t651-yy97.
Genomes to Fields. 2022. Genomes to Fields 2020 dataset. CyVerse Data Commons. doi:10.25739/hzzs-a865.
Gilmour  AR, Cullis  BR, Verbyla  AP, Verbyla  AP. 1997. Accounting for natural and extraneous variation in the analysis of field experiments. J Agric Biol Environ Stat. 2 (3 ):269.doi:10.2307/1400446.
Gilmour  AR, Gogel  BJ, Cullis  BR, Welham  SJ, Thompson  R. 2015. ASReml User Guide Release 4.1 Functional Specification. Hemel Hempstead, UK: VSN International Ltd.
Gitelson  AA, Kaufman  YJ, Stark  R, Rundquist  D. 2002. Novel algorithms for remote estimation of vegetation fraction. Remote Sens Environ. 80 (1 ):76–87. doi:10.1016/S0034-4257(01)00289-9.
Griffing  B . 1956. Concept of general and specific combining ability in relation to diallel crossing systems. Aust J Biol Sci. 9 (4 ):463–493. doi:10.1071/BI9560463.
Heil  K, Schmidhalter  U. 2017. The application of EM38: determination of soil parameters, selection of soil sampling points and use in agriculture and archaeology. Sensors. 17 (11 ):2540. doi:10.3390/s17112540.29113048
Hoefler  R, González-Barrios  P, Bhatta  M, Nunes  JAR, Berro  I, Nalin  RS, Borges  A, Covarrubias  E, Diaz-Garcia  L, Quincke  M, et al  2020. Do spatial designs outperform classic experimental designs?  J Agric Biol Environ Stat. 25 (4 ):523–552.doi:10.1007/s13253-020-00406-2.
Hunt  ER, Doraiswamy  PC, McMurtrey  JE, Daughtry  CST, Perry  EM, Akhmedov  B. 2013. A visible band index for remote sensing leaf chlorophyll content at the canopy scale. Int J Appl Earth Obs Geoinf. 21 :103–112. doi:10.1016/j.jag.2012.07.020.
Kirkpatrick  M, Lofsvold  D, Bulmer  M. 1990. Analysis of the inheritance, selection and evolution of growth trajectories. Genetics. 124 (4 ):979–993. doi:10.1093/genetics/124.4.979.2323560
Lima  DC, Aviles  AC, Alpers  RT, Perkins  A, Schoemaker  DL, Costa  M, Michel  KJ, Kaeppler  S, Ertl  D, Romay  MC, et al  2023. 2020–2021 field seasons of maize GxE project within the genomes to fields initiative. BMC Res Notes. 16 (1 ):219. doi:10.1186/s13104-023-06430-y.37710302
McFarland  BA, AlKhalifah  N, Bohn  M, Bubert  J, Buckler  ES, Ciampitti  I, Edwards  J, Ertl  D, Gage  JL, Falcon  CM, et al  2020. Maize genomes to fields (G2F): 2014–2017 field seasons: genotype, phenotype, climatic, soil, and inbred ear image datasets. BMC Res Notes. 13 (1 ):71. doi:10.1186/s13104-020-4922-8.32051026
Meuwissen  THE, Hayes  BJ, Goddard  ME. 2001. Prediction of total genetic value using genome-wide dense marker maps. Genetics. 157 (4 ):1819–1829.doi:10.1093/genetics/157.4.1819.11290733
Michel  KJ, Lima  DC, Hundley  H, Singan  V, Yoshinaga  Y, Daum  C, Barry  K, Broman  KW, Robin Buell  C, de Leon  N, et al  2022. Genetic mapping and prediction of flowering time and plant height in a maize stiff stalk MAGIC population. Genetics. 221 (2 ):iyac063. doi:10.1093/genetics/iyac063.35441688
Misztal  I, Tsuruta  S, Strabel  T, Auvray  B, Druet  T, Lee  DH. 2002. BLUPF90 and Related Programs (BGF90). Proceedings of the 7th World Congress on Genetics Applied to Livestock Production, 33: 743–44.
Morales  N, Bauchet  GJ, Tantikanjana  T, Powell  AF, Ellerbrock  BJ, Tecle  IY, Mueller  LA. 2020. High density genotype storage for plant breeding in the chado schema of breedbase. PLoS One. 15 (11 ):e0240059. doi:10.1371/journal.pone.0240059.33175872
Morales  N, Kaczmar  NS, Santantonio  N, Gore  MA, Mueller  LA, Robbins  KR. 2020. ImageBreed: open-access plant breeding web–database for image-based phenotyping. Plant Phenome J. 3 (1 ):e20004. doi:10.1002/ppj2.20004.
Morales  N, Ogbonna  AC, Ellerbrock  BJ, Bauchet  GJ, Tantikanjana  T, Tecle  IY, Powell  AF, Lyon  D, Menda  N, Simoes  CC, et al  2022. Breedbase: a digital ecosystem for modern plant breeding. G3 (Bethesda). 12 (7 ):jkac078. doi:10.1093/g3journal/jkac078.35385099
Patrignani  A, Ochsner  TE. 2015. Canopeo: a powerful new tool for measuring fractional green canopy cover. Agron J. 107 (6 ):2312–2320.doi:10.2134/agronj15.0150.
Pebesma  EJ . 2004. Multivariable geostatistics in S: the Gstat package. Comput Geosci. 30 (7 ):683–691. doi:10.1016/j.cageo.2004.03.012.
Pérez-Valencia  DM, Rodríguez-Álvarez  MX, Boer  MP, Kronenberg  L, Hund  A, Cabrera-Bosquet  L, Millet  EJ, van Eeuwijk  FA. 2022. A two-stage approach for the spatio-temporal analysis of high-throughput phenotyping data. Sci Rep. 12 (1 ):3177. doi:10.1038/s41598-022-06935-9.35210494
Piepho  HP, Möhring  J, Williams  ER. 2013. Why randomize agricultural experiments?  J Agron Crop Sci. 199 (5 ):374–383. doi:10.1111/jac.12026.
Pix4D, S. A.  2017. Pix4Dmapper 4.1 User Manual. Lausanne, Switzerland: Pix4D SA.
Poland  JA, Nelson  RJ. 2011. In the eye of the beholder: the effect of rater variability and different rating scales on QTL mapping. Phytopathology. 101 (2 ):290–298. doi:10.1094/PHYTO-03-10-0087.20955083
Reynolds  D, Baret  F, Welcker  C, Bostrom  A, Ball  J, Cellini  F, Lorence  A, Chawade  A, Khafif  M, Noshita  K, et al  2019. What is cost-efficient phenotyping? Optimizing costs for different scenarios. Plant Sci. 282 :14–22. doi:10.1016/j.plantsci.2018.06.015.31003607
Robbins  KR, Backlund  JE, Schnelle  KD. 2012. Spatial corrections of unreplicated trials using a two-dimensional spline. Crop Sci. 52 (3 ):1138–1144. doi:10.2135/cropsci2011.08.0417.
Rodríguez-Álvarez  MX, Boer  MP, van Eeuwijk  FA, Eilers  PHC. 2018. Correcting for spatial heterogeneity in plant breeding experiments with P-splines. Spat Stat. 23 :52–71. doi:10.1016/j.spasta.2017.10.003.
Rutkoski  J, Benson  J, Jia  Y, Brown-Guedira  G, Jannink  J-L, Sorrells  M. 2012. Evaluation of genomic prediction methods for fusarium head blight resistance in wheat. Plant Genome. 5 (2 ):51–61. doi:10.3835/plantgenome2012.02.0001.
Rutkoski  J, Poland  J, Mondal  S, Autrique  E, Pérez  LG, Crossa  J, Reynolds  M, Singh  R. 2016. Canopy temperature and vegetation indices from high-throughput phenotyping improve accuracy of pedigree and genomic selection for grain yield in wheat. G3 (Bethesda). 6 (9 ):2799–2808. doi:10.1534/g3.116.032888.27402362
Sagan  V, Maimaitijiang  M, Sidike  P, Eblimit  K. 2019. UAV-based high resolution thermal imaging for vegetation monitoring, and plant phenotyping using ICI 8640 P, FLIR Vue Pro R 640, and thermomap cameras. Remote Sens. 11 (3 ):330.doi:10.3390/rs11030330.
Sakurai  K, Toda  Y, Kajiya-Kanegae  H, Ohmori  Y, Yamasaki  Y, Takahashi  H, Takanashi  H, Tsuda  M, Tsujimoto  H, Kaga  A, et al  2021. Time-series multi-spectral imaging in soybean for improving biomass and genomic prediction accuracy. Plant Genome. 15 (4 ):e20244. doi:10.1002/tpg2.20244.
Schaeffer  LR . 2004. Application of random regression models in animal breeding. Livestock Prod Sci. 86 (1–3 ):35–45. doi:10.1016/S0301-6226(03)00151-9.
Selby  P, Abbeloos  R, Backlund  JE, Basterrechea Salido  M, Bauchet  G, Benites-Alfaro  OE, Birkett  C, Calaminos  VC, Carceller  P, Cornut  G, et al  2019. BrAPI—an application programming interface for plant breeding applications. Bioinformatics. 35 (20 ):4147–4155. doi:10.1093/bioinformatics/btz190.30903186
Smith  AB, Cullis  BR, Thompson  R. 2005. The analysis of crop cultivar breeding and evaluation trials: an overview of current mixed model approaches. J Agric Sci. 143 (6 ):449–462. doi:10.1017/S0021859605005587.
Sun  D, Robbins  K, Morales  N, Shu  Q, Cen  H. 2021. Advances in optical phenotyping of cereal crops. Trends Plant Sci. 27 (2 ):191–208. doi:10.1016/j.tplants.2021.07.015.34417079
Szeg  G . 1975. Orthogonal Polynomials. 4th edition. Vol. XXIII . Rhode Island: American Mathematical Society Colloquium Publications.
Taghavi Namin  S, Esmaeilzadeh  M, Najafi  M, Brown  TB, Borevitz  JO. 2018. Deep phenotyping: deep learning for temporal phenotype/genotype classification. Plant Methods. 14 (1 ):66. doi:10.1186/s13007-018-0333-4.30087695
Thorp  KR, Thompson  AL, Harders  SJ, French  AN, Ward  RW. 2018. High-throughput phenotyping of crop water use efficiency via multispectral drone imagery and a daily soil water balance model. Remote Sens. 10 (11 ):1682. doi:10.3390/rs10111682.
Van der Werf  JHJ, Goddard  ME, Meyer  K. 1998. The use of covariance functions and random regressions for genetic evaluation of milk production based on test day records. J Dairy Sci. 81 (12 ):3300–3308. doi:10.3168/jds.S0022-0302(98)75895-3.9891276
van Eeuwijk  FA, Bustos-Korts  D, Millet  EJ, Boer  MP, Kruijer  W, Thompson  A, Malosetti  M, Iwata  H, Quiroz  R, Kuppe  C, et al  2019. Modelling strategies for assessing and increasing the effectiveness of new phenotyping techniques in plant breeding. Plant Sci. 282 :23–39. doi:10.1016/j.plantsci.2018.06.018.31003609
Van Es  HM, Van Es  CL. 1993. Spatial nature of randomization and its effect on the outcome of field experiments. Agron J. 85 (2 ):420–428. doi:10.2134/agronj1993.00021962008500020046x.
VanRaden  PM . 2008. Efficient methods to compute genomic predictions. J Dairy Sci. 91 (11 ):4414–4423. doi:10.3168/jds.2007-0980.18946147
White  JW, Andrade-Sanchez  P, Gore  MA, Bronson  KF, Coffelt  TA, Conley  MM, Feldmann  KA, French  AN, Heun  JT, Hunsaker  DJ, et al  2012. Field-based phenomics for plant genetics research. Field Crops Res. 133 :101–112. doi:10.1016/j.fcr.2012.04.003.
Wiesner-Hanks  T, Wu  H, Stewart  E, DeChant  C, Kaczmar  N, Lipson  H, Gore  MA, Nelson  RJ. 2019. Millimeter-level plant disease detection from aerial photographs via deep learning and crowdsourced data. Front Plant Sci. 10 :1550. doi:10.3389/fpls.2019.01550.31921228
Wilkinson  MD, Dumontier  M, Jan Aalbersberg  IJ, Appleton  G, Axton  M, Baak  A, Blomberg  N, Boiten  J-W, da Silva Santos  LB, Bourne  PE, et al  2016. The FAIR guiding principles for scientific data management and stewardship. Sci Data. 3 (1 ):160018. doi:10.1038/sdata.2016.18.26978244
Xu  Y . 2016. Envirotyping for deciphering environmental impacts on crop plants. Theor Appl Genet. 129 (4 ):653–673. doi:10.1007/s00122-016-2691-5.26932121
