Genetic Architecture of Phenomic-Enabled Canopy Coverage in Glycine max.

Digital imagery can help to quantify seasonal changes in desirable crop phenotypes that can be treated as quantitative traits. Because limitations in precise and functional phenotyping restrain genetic improvement in the postgenomic era, imagery-based phenomics could become the next breakthrough to accelerate genetic gains in field crops. Whereas many phenomic studies focus on exploratory analysis of spectral data without obvious interpretative value, we used field images to directly measure soybean canopy development from phenological stage V2 to R5. Over 3 years, we collected imagery using ground and aerial platforms of a large and diverse nested association panel comprising 5555 lines. Genome-wide association analysis of canopy coverage across sampling dates detected a large quantitative trait locus (QTL) on soybean (Glycine max, L. Merr.) chromosome 19. This QTL provided an increase in yield of 47.3 kg ha-1 Variance component analysis indicated that a parameter, described as average canopy coverage, is a highly heritable trait (h2 = 0.77) with a promising genetic correlation with grain yield (0.87), enabling indirect selection of yield via canopy development parameters. Our findings indicate that fast canopy coverage is an early season trait that is inexpensive to measure and has great potential for application in breeding programs focused on yield improvement. We recommend using the average canopy coverage in multiple trait schemes, especially for the early stages of the breeding pipeline (including progeny rows and preliminary yield trials), in which the large number of field plots makes collection of grain yield data challenging.

INFORMATION on spectral reflectance obtained from digital imagery can indicate a plant’s chemical composition and physical properties. Imagery also provides useful information about plant architecture (Yol ). Image data can produce multiple time series and spectral indices (White ), for which limiting factors include the frequency with which pictures are taken during the growing season, and the capabilities of the camera or sensor in a given wavelength spectrum (Andrade-Sanchez ). Recent phenomic studies employing imaging have focused on multi- and hyper-spectral capabilities (Yang ; Kumar ), with minimal investigation of time series of images for capturing crop development and canopy architecture (Schunk ). Yet plant architecture comprises a set of feature-driven traits that can contribute to genetic improvement of field crops. For instance, dwarfing (an architectural trait) is well-known for contributing to the yield increases attained from the Green Revolution (Hammer ; Swaminathan 2014).

Imagery-based field phenotyping with high-throughput platforms has the potential to drive the discovery of novel traits, and to facilitate routine quantification of architectural, developmental, physiological, and phenological crop characteristics (Araus and Cairns 2014). Because it allows inexpensive evaluation of large numbers of field plots with multiple measures in a short period of time (White ), such information provides an efficient framework to tackle phenotyping bottlenecks in plant breeding (Furbank and Tester 2011; Araus and Cairns 2014; Kumar ).

In contrast to point- or plant-based spectral reflectance, image analysis allows direct measurement of spatial or field-based traits that are known to be valuable, such as canopy coverage and vegetation indices, which are possible to collect from the field with high-throughput platforms, such as unmanned aircraft systems (UAS) (Cabrera-Bosquet ; Liebisch ). Additionally, acquisition of imagery can be a robust approach, regardless of its collection platform, that is less sensitive to technical challenges and prevailing environmental conditions than other data acquisition methods. Thus, characterization and selection of phenotypes acquired from image analysis represents a low-hanging fruit for genetic improvement of yield potential in crops.

Increased crop yields often associate with greater biomass production (expressed as the product of total incident radiation throughout the duration of crop growth), the amount of light interception (LI) by crop canopies, and the conversion of chemical energy into plant dry matter (Richards 2000). Improved LI from rapid canopy development is a significant contributor to increased biomass, total photosynthesis, and yield, which makes it an important target when selecting for genetically superior crops. Richards (2000) proposed that genetic selection for improved LI has been occurring inadvertently since the beginning of crop domestication through the selection of young, competitive seedlings with fast crop growth rates, with this resulting in rapid canopy establishment, and, therefore, greater LI during early vegetative growth. Many canopy parameters play dominant roles in determining source-sink ratios (Board and Harville 1993) along with light-use efficiency of the plant’s photosynthetic machinery (Tharakan ).

From 1924 to 2012, on-farm soybean yields increased ∼23 kg ha−1 annually, due to both improved genetics and improved agronomic practices (Rincker ; Specht et al. 2014; Suhre ). Since the introduction of genomic tools in the past two decades, the rate of yield increase has grown to 66 kg ha−1 per year, where 50–66% of gain is associated with genetic improvement (Board and Kahlon 2011; Rowntree ). Recently, evaluation of soybean cultivars released between 1923 and 2007 determined that the greater biomass production and yield potential of newer cultivars is driven by linear increases in LI, efficiency of radiation use, and partitioning efficiency (Koester ). At present, much of the existing genetic variation in photosynthetic efficiency and related processes is confined to undomesticated germplasm (Sherman-Broyles ).

Purcell (2000) developed a precise and rapid technique using ground-based digital imagery for evaluating seasonal variation in soybean canopy coverage. The methodology is based on the binary classification of individual pixels from plot images into canopy or background, thereby providing a direct quantification for the percent of canopy coverage on a continuous scale (Purcell 2000; Karcher and Richardson 2005). This research established a one-to-one relationship between fractional canopy coverage (determined through image analysis) and canopy light interception measurements acquired using a line quantum sensor, thereby facilitating the use of digital imagery to evaluate seasonal canopy interception of photosynthetically active radiation (PAR) (Purcell ; Edwards ). While these findings greatly improved the efficiency with which we can assess LI, this method is impractical for use in large experiments that require phenotyping thousands of genetically distinct lines to determine gene–phenotype associations for complex traits (Myles ).

Rapid canopy development in soybean optimizes growth dynamic parameters, providing a foundation for a greater biomass accumulation during the season, and ultimately greater grain yield potential (Hall 2015); it also favors early-season weed suppression (Jannink ; Fickett ). Canopy coverage can be observed for the entire season, and it is likely that observed variation has a genetic basis, and could be a target for genetic improvement. The objectives of this study were (1) to characterize imagery-based early-season canopy coverage in soybean as a trait with quantitative genetic control, (2) to describe the genetic architecture of canopy coverage via genome-wide association, and (3) to provide insight into the application of UAS quantification of canopy parameters for applications in soybean breeding aimed at improving yield.

The application of marker-assisted selection (MAS) in plant breeding represented a critical step in selection practices driven by realized genetic information (Muir 2007; Hayes ). The applications of genomic tools evolved to whole-genome selection techniques, effective in selecting the favorable alleles in small-effect quantitative trait loci (QTL) that govern low-heritability traits (Heffner ). However, poor phenotypic data quality is a major limitation for effective use of genomic prediction models (Cobb ; de los Campos ). More precise phenotypic information is necessary to accelerate the annual rate of genetic gain in field crops in the postgenomic era (Henryon ). Nevertheless, phenomic platforms can be successful in crop improvement only if they satisfy basic requirements that include high-throughput at a low cost (White ; Basu et al. 2015), and the ability to provide meaningful, heritable, and measurable trait information (Cobb ; Liebisch ).

Imagery acquired from ground and aerial platforms was used to characterize early-season canopy development for 3 years in a next-generation soybean genetic population, for which we also assessed agronomic performance. We developed methodology to automatically extract canopy data from the best quality aerial image of each plot, rather than using mosaics. Percentage canopy pixels was quantified from over 70,000 ground images and over 2000 aerial images. Methods necessary for the analyses include spatial statistics, generalized mixed linear models, and logistic growth function to model longitudinal data of canopy coverage for each genetic line. Genomic data for the population helped determine the genetic architecture of the newly-described soybean trait. Quantitative genetic properties were evaluated to assess the application of high-throughput field-based canopy phenotyping to increase the rate of genetic gain in soybean breeding pipelines.

Materials and Methods

Germplasm

The SoyNAM population (www.soynam.org) is a nested association mapping panel consisting of ∼5600 F5-derived recombinant inbred lines (RILs) that originated as sets of 140 RILs descended from each of 40 biparental matings of 40 founder parents with one common, high yielding, hub parent (IA3023). The RILs exhibit maturities ranging from late MG II to early MG IV. The 40 founder parents include 17 cultivars or elite public breeding lines contributed by soybean breeders across the North Central Region, 15 unique breeding lines developed to have diverse ancestry, and eight diverse plant introductions selected for their performance in water-limited environments. Visual canopy field images of the SoyNAM parents are available at http://www.soybase.org/SoyNAM/imagebrowser.php. Genotyping of the lines employed an Illumina SoyNAM BeadChip SNP array designed for this population, using 5305 single nucleotide polymorphism (SNP) markers identified from the genomic sequences of all 41 parental lines. We imputed missing SNP locus calls using random forest, and removed SNPs with a minor allele frequency <0.15 (Xavier ) using the R package NAM (Xavier ). This left a final total of 4077 SNPs for the association analysis.

Field design

In 2013, 2014, and 2015, we collected phenotypic data from the SoyNAM population in West Lafayette, Indiana. From 2012 to 2014, the experiment used a modified augmented design, and an augmented complete block design in 2015. Line planting took place on May 20, 2013; May 24, 2014; and May 23, 2015, at the Purdue University Agronomy Center for Research and Education (40°28′20.5′′N 86°59′32.3′′W). Experimental units were based on two-row plots (0.76 × 2.90 m) at a density of ∼35 plants m−2. The 2012–2014 experiments grew all the SoyNAM entries, but the 2015 experiment grew only the six families (NAM5, NAM 9, NAM12, NAM15, NAM24, and NAM40) with the largest phenotypic variance for the traits of interest. The 2015 experiment experienced mild flood damage.

Canopy coverage

The metric we used to measure canopy coverage was the percentage of image pixels that were classified as canopy pixels (Figure 1). We took images at regular intervals from 2 to 8 weeks after planting, which spanned soybean phenological stages V2–R5 (Fehr ). The timeframe of data collection coincides with all but the final R5–R7 dry matter accumulation period (Egli and Leggett 1973). Collection of image data used ground-based (2013–2014) and aerial (2014–2015) platforms and images underwent a classification procedure to determine the number of pixels showing the soybean canopy.

Figure 1

Example canopy imagery of a single plot, used to calculate a percentage canopy coverage on a given sampling date. (A, B) From aerial (above; A) or ground (below; B) platforms, with raw (left) and classified (right) imagery.

Classification of ground-based images used SigmaScan Pro software according to the methods of Karcher and Richardson (2005). The Appendix (Supplemental Material, Figure S1) provides examples of classified ground imagery across the seasons. Classification of aerial images used a binomial model implemented in ENVI 5.0 and calibrated with manually labeled training data.

Using a logit link function for each field plot, we fit canopy coverage observations from multiple sampling dates into a logistic model, with the canopy coverage (cc) for any point in time (t) represented aswhere β0 and β1 are the logistic regression coefficients for any given field plot. This function provided daily canopy coverage projections that were used to estimate canopy development parameters. The single-trait representation of canopy in this study is the average canopy coverage (ACC), defined aswhere N is the number of observed days, and t0 and tn represent the first and last day of canopy coverage data collection, respectively. In essence, the ACC value is an arithmetical mean of multiple seasonally observed values of canopy coverage.

Ground-based imagery

Collection of images used to infer canopy coverage employed the methodology described by Purcell (2000). We took Red-Green-Blue (RGB) pictures (680 × 480 pixels) at the center of each plot, ∼1.5 m above the ground, at a 30° angle from the horizon, using a pocket camera (Canon PowerShot A400IS), then cropped these to 560 × 480 pixels to avoid capturing canopy from neighboring plots.

Aerial imagery

In 2014, we collected aerial images using a Precision Hawk UAS equipped with a Nikon 1-J1 digital camera (3872 × 2592 pixels). It collected individual images at an altitude of 100 m, yielding a spatial resolution of 3.25 cm pixel−1. In 2015, it collected images at 50 m using a Nikon 1-J3 digital camera (4608 × 3072 pixels), yielding a spatial resolution of 1.50 cm pixel−1. These images had a 70% lateral and forward overlap, which allowed image stitching using Pix4Dmapper software. Image stitching outputs, including estimated camera positions and orientations, and plot map coordinates, helped to extract ortho-rectified image mosaics, digital surface models, camera positions, and orientations for each aerial photo, and camera model parameters. These outputs were registered to map coordinates (easting, northing, and altitude) using ground control points. We obtained the map coordinates of the individual plots by gridding the image mosaics. The map coordinates of the plots, positions, and orientations of the camera, camera model parameters, and a collinearity relationship, were used to automatically identify which plots appeared in which images and extract ortho-rectified images of the plots directly from the aerial images (Kasser 2002).

Correspondence between platforms

Ganey and Block (1994) reported that canopy coverage estimates from various methods are directly comparable. In 2014, we collected canopy data from both the ground and the UAS platforms, and calculated ACC values for both platforms. A paired data set comparison provided evidence of suitably sufficient correspondence between the two platforms (Figure 2), so we elected to collect data from the UAS platform only in 2015.

Figure 2

Correspondence between the average canopy coverage measured from ground or aerial imagery from 2014.

Statistical model

Correlation and genome-wide association analysis used genetic values obtained from best linear unbiased prediction estimators, rather than observed phenotypes of canopy coverage from individual time points. The following mixed linear model was used to provide the genetic valueswhere y is the vector of observed phenotypes, μ is the intercept, and f(x) is a nonlinear function that accounts for the spatial heterogeneity of field variation within the blocks imposed by the field design, where f(x) is computed as the average phenotypic value of neighbor plots (Lado ). Z is the incidence matrix of environment, u is the vector of regression coefficients of environment effects, W is the incidence matrix of genotypes, g is the vector of genetic values, and e is the vector of residuals. Environment and genotypes were treated as random effects. Random coefficients were assumed to be normally and independently distributed as and

Association analysis

Genome-wide association analysis used the random effect model designed for multi-parental populations (Wei and Xu 2015), implemented in the function gwas2 from the R package NAM (Xavier ). Analysis were performed for individual canopy coverage measurement days spanning 14–56 days after planting, using a random linear effect model:where g is the vector of genetic values of canopy coverage for a given point in time fitted in Equation 3, μ is the intercept, X is the incidence matrix of alleles generated from the marker data and family information, α is the vector of regression coefficients corresponding to the allele effects, ψ corresponds to the polygenic effect that accounts for population structure, and ε is the vector of residuals. Random effect coefficients were assumed to be normal as and where K represents the genomic relationship matrix.

We evaluated statistical significance using a likelihood ratio test by comparing the restricted log-likelihood of the model with the marker included () to the reduced model (), not including the marker (). Thus,In the random model, the LRT follows a mixture of chi-squared and binomial distributions (Xavier ), with P-values computed using a chi-squared distribution with 0.5 degrees of freedom. The Bonferroni threshold that accounts for false positives under multiple hypothesis testing () was used to define which markers were associated with the trait of interest. The expected threshold for 4077 SNPs markers was estimated 4.91 –log (P-value).

Phenotyping of agronomic traits

Grain yield was measured as grams of seed per plot, converted to kg ha−1 using harvest-timed seed moisture to adjust all plot values to 13% seed moisture. We collected maturity data twice a week as the number days after planting (DAP), back and forward scoring plots that flowered and matured between the intervals. The criterion for a plot to achieve maturity (R8) was 50% of the plants having 95% of their pods mature (Fehr ).

Variance component analysis

We computed Pearson and Spearman correlations from the phenotypic data using built-in R functions. A multivariate mixed linear model estimated covariance components for each of the three traits: grain yield, days to maturity, and average canopy coverage. The linear model was presented on Equation 3, but environment was treated as fixed effect. Model fitting used the AI-REML algorithm (Gilmour ) implemented in AIREMLF90 (Misztal ) with the covariance structure set up as follows:Here K represents the realized genomic relationship matrix. Matrices and represent the genetic and residual covariance matrices among the three traits. These covariance components then helped to estimate the heritability of average canopy coverage, grain yield, and days to maturity, as well as the genetic correlations among these traits.

Data availability

Phenotypes and genotypic data are available in the R package NAM (https://CRAN.R-project.org/package=NAM). Load the data using the following command: data(met, package = “NAM”).

Results and Discussion

Canopy behavior

Individual lines showed the potential to attain full canopy coverage (100%) as early as 50 DAP, while others did not reach full coverage within 64 DAP (Figure 3). In the latter case, lines seemed incapable of closing the canopy entirely in a production system employing a row spacing of 76 cm. Reports indicate that yield and canopy traits are under similar genetic control when grown using narrower row spacing, with slight changes in heritability (Weaver and Wilcox 1982; De Bruin and Pedersen 2008).

Figure 3

Development of canopy coverage in the SoyNAM population described as logistic curves of canopy development for individual soybean lines over two or three seasons. Each curve describes an individual soybean line. The white dotted curve represents the mean curve.

Genome-wide association

Six genomic regions, on soybean chromosomes 1, 5, 6, 9, 10, and 19, associated with canopy coverage when measured or estimated across all days (Figure 4b). The largest effect QTL detected was on chromosome 19, and it associated with canopy coverage on all days within the survey period (Figure 4a). The second largest effect QTL was located on chromosome 10, and was significant for approximately two-thirds of the survey period, from 14 to 26 DAP. Table 1 summarizes the QTL detected, the period during which these QTL were significant, and their average allelic effects on ACC, grain yield, and maturity.

Figure 4

Results of the genome-wide association study for canopy coverage. Genomic regions significantly associated with early-season soybean canopy coverage for 2 or 3 years observed or estimated for (A) each day from 14 to 56 days after planting, and (B) across all days.

Table 1

SNPs significantly associated with early-season soybean canopy coverage and the number of DAP during which they were significant

SNP	Period of Significant Association (DAP)	GY (kg ha−1)	R8 (DAP)	ACC (%)
Gm01_50911939_C_T	56	−26.04	−1.12	−0.36
Gm05_37467797_A_G	53–56	−44.71	0.08	−0.65
Gm06_14104090_T_C	56	99.58	0.50	0.55
Gm09_4034850_C_T	14–35, 59–64	12.24	−0.60	−0.44
Gm010_44120764_T_C	14–23	5.95	−0.59	0.04
Gm010_44630777_C_A	14–26	51.61	1.06	0.33
Gm019_1586092_T_C	14–64	47.30	−0.24	1.34

Mean change in phenotypic values associated with each QTL in terms of grain yield (GY) expressed as kilograms per hectare, number of DAP to maturity (R8), and the percentage of ACC, calculated as the allelic substitution of the homozygous founder parent genotype for the homozygous IA3023 genotype.

The existence of large-effect QTL along with several nearly significant associations indicate that ACC is controlled by major genes in the soybean genome. Another canopy trait, early-season canopy height, was reported to be heritable ( = 0.64) and genetically correlated to weed-suppression ability () (Jannink ).

Genome-wide association presented by day (Figure 4a) shows that the chromosome 19 QTL was significant throughout the survey period. Other associations were significant either earlier or later in the season, or when the canopy was nearly closed (Table 1). Early- and late-season canopy growth are controlled through independent genetic mechanisms and can be dissected into distinct traits (Jannink ).

Assuming pleiotropy, the estimated increase in grain yield associated with the positive allele of the chromosome 19 QTL was 47.30 kg ha−1. Importantly, this increase did not also increase days to maturity (−0.24 days).

Quantitative genetic properties of average canopy closure

Variance components provide an insight about the level of genetic and environmental control of the traits of interest, necessary to estimate heritability (Table 2). At multivariate level, (co)variances indicate the nature of the relationship between multiple traits, uncovering whether traits are more associated at genetic or environmental level (Table 3).

Table 2

Variance components of the traits considered in this study

	GY	R8	ACC
Var(G)	78.80	208.36	25.54
Var(E)	56.12	9.39	8.10
h^2	0.58	0.96	0.76

Genetic (G) and environmental (E) variances and heritabilities for soybean grain yield (GY), and days to maturity (R8), and ACC.

Table 3

Phenotypic Pearson (P), phenotypic Spearman (S), genetic (G), and environmental (E) correlations among GY, ACC and days to maturity (R8)

	GY-ACC	R8-ACC	GY-R8
Cor(P)	0.63	0.31	0.42
Cor(S)	0.70	0.38	0.46
Cor(G)	0.88	0.77	0.72
Cor(E)	0.18	−0.06	0.23

For ACC to be an effective predictor of grain yield potential in the early stages of breeding pipelines, the phenotypic correlations, and, more importantly, the genetic correlations, must be high between the two traits. The genetic correlation between grain yield and ACC is quite high (Table 3), but so are the correlations between days to maturity and ACC, and between grain yield and days to maturity. However, as noted above, the main allelic effect of the chromosome 19 QTL affects ACC and yield, without negatively impacting days to maturity.

Selection for canopy coverage

A faster rate of canopy coverage development (achieved through manipulation of crop management practices) increases light interception, suppresses weeds, and has other beneficial properties (Purcell 2000; Jannink ; Martius ; Campillo ). As a target for artificial selection, the trait presents attractive genetic properties, such as high heritability and strong genetic correlation with grain yield that allows for indirect selection or prediction of grain yield potential.

For indirect selection to be effective, the response of grain yield potential when selecting for ACC must be greater than the response of selecting for yield directly. Very few cases of successful indirect selection have been reported (Bernardo 2010). The indirect selection equation indicates that efficient indirect selection for yield at the same selection intensity is possible with ACC because the correlated response (CR) to selection for grain yield is 14% greater than the direct response (R):In the equation above, represents the heritability of the secondary trait (i.e., ACC), is the genetic correlation between primary and secondary traits, and is the heritability of the primary trait (i.e., grain yield potential). Thus, assuming identical selection intensities, indirect selection of grain yield using average canopy coverage would be expected to be better than direct selection for yield. Although the canopy traits reported by Jannink did not show maturity-mediated tradeoffs, use of selection indices including both ACC and days to maturity may be necessary to mitigate indirect increases in maturity.

The use of average canopy coverage to predict lines with high grain-yield potential when collection of accurate yield data are not feasible or is challenging. This is often the case during early generations in a breeding pipeline when yield data are of low quality due to low seed number, or when limited resources are available for the evaluation of many lines, as in preliminary yield trials. ACC may also be useful in selection indices that include grain yield and other valuable traits.

UAS phenotyping and phenomics

From the perspective of genetic improvement of crops, traits with known properties that can be collected in situ using UAS platforms are a low-hanging fruit, in particular canopy characteristics and vegetation indices (Cabrera-Bosquet ; Liebisch ). In addition to being robust, simple, and cost-effective, these traits are likely more valuable than measurements made in controlled environments for physiological characteristics that cannot be translated into field phenotypes, and which seldom provide economic advantages (White ; Von Mogel 2013). Here, we have shown canopy coverage to be a highly heritable trait associated with grain yield potential, possibly the most promising phenomic trait reported so far.

Accounting for environmental noise

This study used spatial statistics in its model, which provided a covariate term that accounted for micro-environmental field variation, or environmental noise, on a plot-by-plot basis (Lado ). Without this adjustment, the heritability of ACC was estimated to be 0.46 (data not shown), though the magnitude of the genetic correlation of ACC with grain yield remained the same. Previous studies (Piepho ; Oakey ) described the importance of considering spatial structure among field observations to achieve proper estimation of statistically-estimated trait values, or breeding values. Many authors report the challenge of differentiating between signal and noise in phenomic data collected from the field (Cobb ; Araus and Cairns 2014; Basu et al. 2015). We conclude that analysis of field phenomic traits should not ignore spatial information.

Breeding values incorporating phenomic data

Research is beginning to explore the optimal application of phenomic traits in genetic improvement of crop plants, and how selection schemes and prediction models should accommodate phenomic information. Phenomic quantifications from field imagery are indicators of plant health and development (Liebisch ), but inclusion of such measurements as covariates in genomic models can be controversial (Valente ) because they share both a genetic and an environmental correlation with agronomic traits (Table 3). Instead, phenomic traits might be more suitably used in multi-trait models, to improve the accuracy of estimated breeding values by exploiting genetic correlations, and, therefore, preserve the genetic nature of the association between the target trait and phenomic trait. We showed that canopy coverage displays a high correlation with valuable soybean traits (Table 3), that favors multivariate models and enhances the accuracy of breeding values.

Supplementary Material

Supplemental material is available online at www.genetics.org/lookup/suppl/doi:10.1534/genetics.116.198713/-/DC1.

Click here for additional data file.