Development and Validation of a Phenotyping Computational Workflow to Predict the Biomass Yield of a Large Perennial Ryegrass Breeding Field Trial

Abstract

Increasing dry matter yield (DMY) is the most important objective in perennial ryegrass breeding programs. Current yield assessment methods like cutting are time-consuming and destructive, non-destructive measures such as scoring yield on single plants by visual inspection may be subjective. These assessments involve multiple measurements and selection procedures across seasons and years to evaluate biomass yield repeatedly. This contributes to the slow process of new cultivar development and commercialisation. This study developed and validated a computational phenotyping workflow for image acquisition, processing and analysis of spaced planted ryegrass and investigated sensor-based DMY yield estimation of individual plants through normalized difference vegetative index (NDVI) and ultrasonic plant height data extraction. The DMY of 48,000 individual plants representing 50 advanced breeding lines and commercial cultivars was accurately estimated at multiple harvests across the growing season. NDVI, plant height and predicted DMY obtained from aerial and ground-based sensors illustrated the variation within and between cultivars across different seasons. Combining NDVI and plant height of individual plants was a robust method to enable high-throughput phenotyping of biomass yield in ryegrass breeding. Similarly, the plot-level model indicated good to high-coefficients of determination (R²) between the predicted and measured DMY across three seasons with R² between 0.19 and 0.81 and root mean square errors (RMSE) values ranging from 0.09 to 0.21 kg/plot. The model was further validated using a combined regression of the three seasons harvests. This study further sets a foundation for the application of sensor technologies combined with genomic studies that lead to greater rates of genetic gain in perennial ryegrass biomass yield.

Introduction

Increasing biomass yield is the most important trait to improve during the breeding of perennial ryegrass (Smith et al., 2001; McDonagh et al., 2016; Herridge et al., 2018; Ghamkhar et al., 2019). However, biomass yield is a complex trait, which varies with the number and density of tillers, regrowth after defoliation and growth habit; across a range of seasons, environments, and plant age (Yates et al., 2019). Furthermore, the early stages of perennial ryegrass breeding programs depend on the assessment of populations based on large numbers of genotypes planted as spaced plants or small sward plots in the field (Lootens et al., 2016; Ghamkhar et al., 2019). These assessments involve multiple measurement and selection procedures across seasons and years to repeatedly evaluate biomass yield (Leddin et al., 2018). Biomass yield assessment methods are time-consuming, mainly relying on visual scores or destructive harvesting, drying and weighing of samples. This contributes to the slow process of new cultivar development and commercialization. The application of sensor-based high-throughput phenotyping (HTP) technologies on aerial and ground-based mobile platforms have the potential to offer a non-destructive, rapid and efficient method to assess biomass yield in the field (Inostroza et al., 2016; Gebremedhin et al., 2019b; Ghamkhar et al., 2019).

The possibility for the application of HTP for automated biomass assessments on an objective and accurate basis has been demonstrated for forage grasses (Walter et al., 2012; Barrett et al., 2015). In other species, various remote sensing methods have been calibrated and validated to conduct measurements at a single plant, plot or experimental units, paddock and landscape scales (Rahman et al., 2014; Hunt et al., 2015; Inostroza et al., 2016; Punalekar et al., 2018). Satellite-based remote sensing methods have been used to assess grassland biomass at a landscape scale (Barrachina et al., 2015). However, low temporal and spatial resolution limits the applicability of satellite remote sensing for breeding (Tattaris et al., 2016). For instance, satellite imaging tools have spatial resolution ranging from 1.24 to 260 m (Chawade et al., 2019), which limits their application on small experimental plots. Proximal sensing platforms have the advantage of capturing phenotypic data from thousands of genotypes at high spatial and temporal resolution (Haghighattalab et al., 2016). The platforms equipped with various sensors may function to capture information of multiple traits from field plots and individual plants.

A variety of aerial and ground-based HTP platforms have been developed for non-destructive biomass estimation of wheat (Yue et al., 2017; Jimenez-Berni et al., 2018), barley (Bendig et al., 2014; Brocks and Bareth, 2018), triticale (Busemeyer et al., 2013), and rice (Devia et al., 2019). Specifically, small grassland experiments (10 sampled plots of 1 m² each), demonstrated non-destructive biomass estimations of timothy and ryegrass from plant height (Rueda-Ayala et al., 2019). Other experiments of biomass estimation on ryegrass (Borra-Serrano et al., 2019) and a mixture of timothy and meadow fescue (Viljanen et al., 2018) were demonstrated using aerial platforms. In these studies, vegetative indices, plant height and photogrammetric features were validated, and the models from these surrogates were used to develop biomass estimation models under plots and small paddock level (n ≤ 96 sampled plots) experiments. Thus, the general application of aerial and ground-based phenotyping methods in these experiments showed there is a potential to increase throughput and hence to improve the precision to individual plants level biomass yield estimation. Therefore, development and validation of a phenotyping computational workflow of data acquisition, processing and analysis could be applied to implement automated biomass yield estimation in large-scale ryegrass breeding program. We have recently shown that normalized difference vegetative index (NDVI) and plant height correlated with perennial ryegrass biomass yield in four seasons (Gebremedhin et al., 2019a). The results showed combining NDVI and plant height to be a robust method to enable high-throughput phenotyping of biomass yield in a ryegrass breeding program.

The objective of the study was to validate a computational phenotyping workflow for image acquisition, processing and analysis of spaced-planted perennial ryegrass to estimate the biomass yield of 48,000 individual plants through NDVI and plant height data extraction. The workflow resulted in the accurate, non-destructive biomass yield estimation of spaced planted perennial ryegrass across growing seasons in a field environment. The outcome of the study contributes to the application of sensor technologies combined with genomic selection for greater rates of genetic gain in forages.

Materials and Methods

Field Experiment and Plant Material

Advanced perennial ryegrass breeding lines and commercial cultivars were established in a field experiment in June 2016 at Agriculture Victoria Research, Hamilton Centre, Victoria, Australia (37.8464°S, 142.0737°E) (Supplementary Figure 1), as part of a project to implement genomic selection (GS) in perennial ryegrass. Fifty advanced perennial ryegrass breeding lines and commercial cultivars were planted in a randomized complete block design with ten replications of each. Each line or commercial cultivar was coded 1–50 to conceal their identities as requested by the owners of the breeding lines and cultivars. Each replication was considered as a plot and contained three rows of 32 spaced plants each (i.e., 96 plants/plot). The experimental unit was, therefore, a plot of 8 × 1.8 m. The spacing between plants was 25 cm and between rows was 60 cm. The GS experiment contained a total of 48,000 individual plants in 10 blocks. The total area of the field experiment was 8,100 m².

Phenotypic Data Collection

For the non-destructive biomass estimation, ground and aerial-based proximally sensed data were collected a day before each harvest. The harvests were allocated to seasons according to the protocol of Leddin et al. (2018). The field experiment was harvested three times in 2018, one in June 2018 called winter 2018, and two in the late spring season of October 2018 and November 2018 called late-spring2018_1 and late-spring2018_2, respectively. No harvests were conducted in autumn or summer due to slow growth and drought. In these three harvests, two sets of data collection scenarios were conducted at the same time.

Firstly, 480 individual plants from three breeding lines were selected and manually harvested at 5 cm above ground level in each of the seasonal harvests. The manually harvested measurements were used to develop and validate dry matter yield (DMY) estimation model from the combinations of NDVI and plant height values from sensors (Gebremedhin et al., 2019a). The models developed from these individual plants were used to estimate biomass of the 48,000 plants for the winter2018, spring2018_1, and late-spring2018_2 seasonal harvests.

Secondly, mechanical harvesting of the 500 plots of the GS field was conducted using a Gianni Ferrari combine harvester (Gianni Ferrari s.r.l., Reggiolo, Italy) to validate the plot level correlation between estimated and measured DMY. In the late-spring2018_1 and late-spring2018_2 measured DMY data wasn’t collected from 100 plots since these plots were allocated for flowering score experiment. Plant height data from a few plots of the late-spring2018_2 was not received due to technical issues of the Phenorover driving at that specific date.

Aerial Images Acquisition and Processing

Multispectral images acquisition was conducted using a RedEdge-M sensor (RedEdge, MicaSense Inc., Seattle, WA, United States) attached to DJI M100 quadcopter (DJI Technology Co., Shenzhen, China). The RedEdge-M capture images simultaneously at five bands including blue (465–485 nm), green (550–570 nm), red (663–673 nm), red edge (712–722 nm), and near-infrared (820–860 nm). It also has GPS and sensor and incident light sensors. The flight mission was planned by Pix4D Capture software.

Aerial images were collected using the unmanned aerial system (UAS) on a weekly basis over the GS trial site, and three sets of data from 2018 (winter2018, late-spring2018_1, and late-spring2018_2) was used for this analysis. Imaging dates were synchronized with each harvest. Flight operations were conducted under bright, sunny weather conditions to minimize noise from environmental variation. The UAS flight altitude was set at 20 m above ground level, and the flight speed was 6 m/s with 75% side and forward overlap of images. At this flight altitude and speed, the spatial resolution of the images was 2 cm/pixel. The same flight path was followed on each date. Multispectral ground calibration targets (Tetracam Inc., Chatsworth, CA, United States) with known reflectance values (6, 11, 22, and 33%) was used for radiometric calibration by developing regression equation between digital numbers (DN) vs. the reflectance percentage.

Computational Workflow for Data Extraction

Before the first flight, nine ground control points (GCPs), (50 × 50 cm checkered linoleum flooring) were distributed across the field experiment to improve the georeferencing accuracy. The GCPs were used for manual assignment of absolute geographic coordinates to the ortho-mosaic images. For image processing, Pix4Dmapper Pro (version 4.3.31 Pix4D, Lausanne, Switzerland, https://pix4d.com) was used to process raw images into ortho-mosaics. The workflow of analyzing geospatial images of spectral signals (Figure 1) is described below.

FIGURE 1

Ortho-Mosaic Processing, Calibration, and NDVI Extraction

Ortho-mosaic processing have five main steps; including image alignment, image geo-referencing, building dense point clouds, ortho-mosaics, and vegetative indices map, were followed for image processing. For detailed processing steps, please refer to following previous studies (Khan et al., 2018; Dobbels and Lorenz, 2019). Ortho-mosaic images were generated and stored as geo-referenced TIFF files. Assessment accuracy of processed ortho-mosaic images was conducted by assessing the quality report from the Pix4D processing. After the ortho-mosaic processing was completed, reflectance calibration and NDVI extraction were conducted using eCognition. For radiometric calibrations, linear correlations were developed between ortho-mosaic images of DN and their known reflectance values for each band. The equation developed was used to convert DN to corrected reflectance values of the ortho-mosaic TIFF file. The corrected ortho-mosaic TIFF file was used to develop the NDVI extraction workflow in eCognition. NDVI was calculated from the red and near-infrared reflectance values (Khan et al., 2018).

We generated a polygon as the boundary of each row of 32 plants as well as a polygon around the boundary of each individual plant with an identification number (plant ID) from the calibrated NDVI images above. The process comprises of four components:

Perspective Transform

Four corner GCPs in the field were used to form a polygon bounding the area of interest. The area contained in this polygon for each NDVI image was subsequently cropped and rotated into a rectangular shape using a perspective transformation (Figure 2), where the plant rows and blocks were essentially orthogonal to the image axes. Identification of plant row polygons.

FIGURE 2

A set of 1,500 polygons, each of which contains a single row of plants, was generated based on each transformed NDVI image (time point). The best polygon for each plant row, which we refer to as a consensus polygon, was subsequently chosen and applied to all downstream analysis later.

Plant rows (Figure 3D) for each transformed NDVI image were identified using projection methods, whereby the sum of the image pixels along the image columns and image rows was used to produce vertical and horizontal projection profiles. These profiles were then processed (Figure 4) to identify the centers of each trough (low-intensity regions) separating either the plant blocks (Figure 3B) or the plant rows (Figure 3D). Briefly: first, the plant blocks were identified by taking the projection of the image pixels values along the x-axis as shown in Figure 3A where the horizontal projection is seen in the profile to the right of the sub-image. The actual processing steps used to identify the troughs between the plant-blocks are laid out in Figure 4, from which it is seen that only one external parameter is required; that is the expected length of the trial plant rows. Next, once each plant block was isolated (Figure 3C) the individual plant rows were then identified using the same process outlined in Figure 4 in each plant block, but where the projection was along the y-axis, and the plant-row length parameter was replaced with the expected plant width/diameter.

FIGURE 3

FIGURE 4

The edges associated with each plant-block and plant-row (Figure 4D, blue vertical lines) were then used to construct an initial estimate of each plant-row polygon/box (Figure 5A). As can be seen in Figure 5A, the initial estimates of the plant-row boxes underestimated the extent of each box in the vertical direction. The amount of underestimation in this case was deliberate and was controlled by the size of the erosion in the profile processing (Figure 4B). This was done to compensate for the lack of alignment in the actual start and end vertical position of each planted plot-row in each block (Figure 5A). However, this meant the trial plants between two vertically aligned plant-boxes were not assigned to a plant-row box. Therefore, to assign each of these plants to a plant-row-box a box growing procedure was developed as shown in Figure 5B.

FIGURE 5

In the box growing procedure, each initial plant-row box (Figure 5A) is linked with the boxes directly above and or below it. Each initial plant-row box is linked to at most two other boxes. The region between two linked boxes can be defined by a rectangle, with the top edge of the link delineated by the bottom edge of the top box, and the bottom edge of the link delineated by the top edge of the lower box. The sides of the rectangular linking region are defined by the minimum and maximum pixel x-positions taken from the intersection of the x-positions defined by the top and bottom boxes. Therefore, the width of a linked region is defined by the maximum x-position minus the minimum x-position plus 1.

Next, for each of the linking regions, a projection profile was constructed from the row-wise sums of the pixels in each of the rows in the linking region. Let such a profile be represented as f(x), where x∈{1,…,n} and n is the number of image pixel-rows between the linked boxes. f(1) represents the sum of pixels along the top edge of the link and f(n) represents the sum of the pixels along the bottom edge of the link. Let t represent a positive threshold value. The threshold t is set to the average background (non-plant) pixel intensity and is scaled during box growing to the number of pixels between the upper and lower edges of the growing boxes. This threshold can be determined from the image or can be supplied by the user. The details of the box growing procedure are outlined in the following pseudo-code:

Choosing Consensus of Plant Row Polygons

As the plant row polygons were detected for each image separately, we chose a best set of polygons that could be applied to all downstream analysis. In order to choose the most stable set of polygons, the variance of the area of the 1,500 polygons for each transformed NDVI image (time point) was calculated and compared across all time points. The plant row polygons identified from the image captured at 31/03/2017 had the least variance in polygon area and were chosen as the consensus.

Individual Center Detection

To identify single plants within the GS field trial, we applied the following 2-step process.

Firstly, a center point of each individual plant was identified. The NDVI ortho-mosaic TIFF image (captured date: 31/03/2017), from which the consensus of plant row polygons was developed, and subsequently processed in Pix4D was imported into eCognition (eCognition developer 9.2, Trimble, Munich, Germany). The template editor algorithm in eCognition was used to allow for the identification of individual plants. A test template was created by selecting 35 samples (Width: 20; Height 20; Context: 4) which were representative of various sizes of perennial ryegrass plants. The correlation between the generated test template and selected samples was 0.774. A small-region of the NDVI layer was randomly selected to validate the generated template and predict individual plants within the area selected. The identified plants within template were then manually verified by the user as being correct, false, or unsure predictions to allow eCognition to refine the prediction model. After the template validation, test parameters values were obtained [ground truth tolerance (pixels): 3 and threshold: 0.6]. The template applied to the test region correctly identified 281 plants, missed 18 plants and incorrectly identified 2 samples as plants. The generated template was then updated based on the test area and the updated template, which had a sample correlation (0.917) to the samples within the selected area. This template matching algorithm was then applied to the trial area. All points predicted outside of the plant-row box was discarded. Manual visual quality control was performed within each plant-row box and further points added or removed to ensure that it was correct representative of the plants within the image were included.

Individual Polygon Detection

Once all plant-row-boxes were constructed, and all plant centers identified, a simple procedure was used to assign each plant-center to its enclosing plant-row-box. Once all plant-centers were assigned to a plant-row box, individual plant-boxes needed to be constructed around each plant center. This was done by (1) using the midpoint between adjacent plants-centers to create a horizontal line connecting the left and right edge within each plant-row-box and (2) the enclosing boxed region containing a plant-center was then identified and used to determine the coordinates for each plant-box. Consequently, each plant-row polygon was split into 32 individual plant polygons (bottom right Figure 1). Each individual plant polygon was also annotated by an identification number (plant ID) and was used to calculate the NDVI statistics for each plant in eCognition.

Ground-Based Plant Height Data Acquisition and Processing

In the ground-based platform, “PhenoRover,” a Polaris ranger side-by-side vehicle was used for ryegrass morphometric data collection. Details of the PhenoRover platform and data processing method was discussed in our previous paper (Gebremedhin et al., 2019a). In brief, the vehicle has six ultrasonic sonars for plant height measurement attached to a welded steel bar at 0.6 m above-ground. For accurate geolocating of each sensor’s data, a single Real-Time-Kinematic Global Navigation Satellite System (RTK-GNSS) unit was placed on the top of the vehicle. Sensors data was stored on a data logger (CR3000, Campbell Scientific, Inc., Logan, UT, United States). For extracting data, we first projected the latitude and longitude coordinates to Universal Transverse Mercator (UTM) coordinate system and calculated the sensor position followed by matching sensor data to individual plants in the quantitative GIS software interface. The same polygons used to extract NDVI were superimposed on top of each geo-referenced data in the quantitative GIS software interface to extract plant height.

Destructive Biomass Measurements

The aboveground biomass was harvested manually at 5 cm height. Destructive biomass harvest of 480 space-plants was collected following the same procedure described previously (Gebremedhin et al., 2019a). In addition, mechanical cutting of all the GS field experiment plots was conducted to measure total fresh biomass weight (kg/plot) at the row or plot level. In this case, biomass was measured at each seasonal harvest from all 500 plots. This data was used to validate the computational phenotyping workflow with manually derived polygons as well as to compare manually collected biomass with predicted biomass using the workflow.

Data Analysis

The variability distribution of phenotyping data was statistically analyzed using R version 3.6.1, R Development Core Team, Vienna, Austria. Non-destructive biomass estimation of individual plants was described using the combination of NDVI and plant height. Spatial distribution variance maps of NDVI, plant height, predicted and measured DMY in three seasons were compared. Histogram of distribution graphs was prepared in Tableau Software (Tableau2019.4, Seattle, Washington, United States). The predicted DMY of spaced-planted individuals was therefore estimated using the NDVI and height combination (NDVIsq_PH) using the following general and seasonal formulas derived from Gebremedhin et al. (2019a):

Universal Equation:

Seasonal Equations:

The statistical regression between the predicted and measured DMY during each harvest was compared at single plant and plot level. Similarly, a model considering all three harvests of biomass estimation was also developed. Coefficient of determination (R²) and root mean square error (RMSE) were used to evaluate and quantify the accuracies of regression models of phenotyping variables. Normality of the data was verified by checking the residual vs. fitted values of plots to ensure the normality of residuals with constant variance.

Results

Comparison of Automated and Manually Derived Delineated Polygons

In our previous validation experiment, we developed a model to predict the DMY of individually drawn polygon by combining NDVI and plant height (Gebremedhin et al., 2019a). Polygons used to extract data were manually, but it is difficult to apply this into 48,000 plants. As a result, a computational phenotyping workflow was applied to extract NDVI and plant height information. To evaluate the robustness of the automatically derived polygons, the predicted DMY of individual plants from automated and manually derived polygons was compared across three different seasons (Figures 6A–C). The results show a strong coefficient of determination between predicted DMY from automated and manual polygons with R² values ranging 0.82-0.93 for the winter2018 and late-spring2018_1 season and RMSE values ranging from 1.81 to 3.83 g/plant. For the late-spring2018_2, harvest a moderate coefficient of determination (R² = 0.44) of DMY prediction was obtained. These results indicate that automatically built polygons can replace manually derived polygons for plant height and NDVI accurate extraction, and biomass prediction for all 48,000 plants in our field experiment.

FIGURE 6

Seasonal Spatial Variation of Biomass and Other Traits

The plant height, NDVI and predicted DMY variables were retrieved using a computational data processing workflow developed in this experiment (Figures 7A–C). NDVI showed seasonal variability with increased values observed in late-spring2018_1 compared to winter2018 and late-spring2018_2. Most of the plots in the late-spring2018_1 indicated NDVI was at saturation point at the time of harvest. Spatial variations of NDVI values were observed in the middle, top and lower part of the field experiment in all seasons. For instance, the center part of the GS trial site showed explicitly lower NDVI values in winter2018 compared to an upper and lower location within parts of the field experiment. Seasonal plant height and predicted biomass variation occurred in the winter2018 followed by peak accumulation in the late-spring2018_1 then decreased into late-spring2018_2.

FIGURE 7

Distribution of NDVI, Plant Height, and Predicted Biomass Across Growing Seasons

The distribution of NDVI, plant height and predicted biomass across three growing seasons, of the entire population of 48,000 individual perennial ryegrass genotypes was estimated (Figures 8A–C). Figure 8A shows the NDVI distribution variation of individual plants across seasons. The winter2018 season had a wider spread of NDVI values at the time of harvest for most genotypes, with some already reaching saturated NDVI values. In the late-spring2018_1, the narrow spread of NDVI values distribution was skewed toward the right indicating that most of the individual plants reached the NDVI saturation point at the time of harvest. In the late-spring2018_2 harvest, relatively normalized distribution of NDVI values was observed signifying some individual plants were harvested earlier than the 2–3 leaf stage.

FIGURE 8

The plant height of spaced-plants showed a normal distribution across all growing seasons (Figure 8B). Taller plants were observed at the late-spring2018_1 harvest followed by the winter2018 and late-spring2018_2 harvests, respectively. Figure 8C shows the predicted DMY distribution variation of individual plants across seasons. The winter2018 and late-spring2018_2 seasons histogram distributions were similar with left-skewed distributions of narrow-spread values at the time of harvest for most genotypes. The predicted DMY distribution of the late-spring2018_1 harvest was observed to be normally distributed, contrary to the winter2018 and late-spring2018_2 harvests.

Comparison of Seasonal Distribution of NDVI, Plant Height, and Predicted Biomass of Cultivars

Figure 9 displays the NDVI histogram distribution variance of 50 cultivars/breeding lines for each of the three harvests. Cultivars showed large relative differences in growth and development between and within a harvest. For example, cultivars designated with a Cultivar ID number 8, 10, 13, 25, 45, 46, 47, and 48 reached saturation and the overdue harvest was observed in winter2018.

FIGURE 9

Figure 10 displays the plant height variation distribution of 50 cultivars/breeding lines across three seasons. Similar to the NDVI histogram distribution, the plant height variation of each cultivar was obvious at the winter2018 harvest. For instance, cultivars designated with Cultivar ID number 8, 10, 13, 25, 45, 46, 47, and 48 showed higher height values than their counterparts. Similarly, plants with Cultivar ID 5, 9 12, 14, 16, 21, 33, 35, 36, and 41 showed a wider distribution at the late-spring2018_1 harvest. The plant height distribution for the late-spring2018_2 harvest did not show any major histogram structural differences among cultivars.

FIGURE 10

Figure 11 displays the variation of individual plant biomass (predicted) among and within 50 perennial ryegrass cultivars/advanced breeding lines across three growing seasons. Similar to the plant height and NDVI, cultivars designated with a Cultivar ID number 8, 10, 13, 25, 45, 46, 47, and 48 showed higher height values. However, the histogram of predicted DMY for each cultivar at the late-spring2018_1 harvest showed more variation compared to the winter2018 and late-spring2018_2 harvests. Moreover, varietal differences in DMY distribution of some cultivars showed a clear difference at the winter2018 harvest.

FIGURE 11

Plot Level Seasonal Biomass Estimation Model

The automated biomass prediction of spaced individual plants was translated to plot-level estimates. The sum of predicted biomass from spaced planted 96 individual plants within the plot was calculated to determine predicted biomass per plot and compared with the mechanically harvested and measured biomass data. Significantly high coefficients of determinations were observed between measured and predicted DMY (R² = 0.19-0.81) at all three harvests (Figures 12A–C). Subsequently, RMSE values ranged from 0.09 to 0.21 kg/plot was observed at these harvests. When combining all harvests into a single regression of estimated and measured biomass (Figure 13), the best coefficients of determination had a slope, R² and RMSE values of 0.93, 0.94, and 0.21 kg/plot, respectively.

FIGURE 12

FIGURE 13

Discussion

In commercial ryegrass breeding, experiments with thousands of individual space-plants may be used for evaluation and selection (Hayes et al., 2013; Lin et al., 2016). Phenotyping of these genotypes requires repeated evaluation across growing seasons and years. This is essential to select, potentially high yielding genotypes to take into future stages of a breeding program. However, phenotyping of large populations, by the available traditional methods of visual scores, manual cutting, bagging and drying is difficult to implement in terms of time and cost. This study developed and validated ground-and aerial-based platforms equipped with advanced sensors for high-throughput phenotyping of ryegrass biomass yield. This resulted in the development of a computational workflow for image acquisition, processing, and analysis to predict biomass yield based on vegetative indices and plant height from 48,000 ryegrass plants. With this method, it was possible to assess traits including predicted biomass of 48,000 individual plants in 6 h from data acquisition to data extraction. This method may provide a robust method for ranking of genotypes and make faster progress in breeding programs.

Manually drawing thousands of polygons has been indicated difficult and considered as time-consuming and automation of polygons for 38,000 plots was expected to increase the throughput and accuracy of traits extraction (Shi et al., 2016). A coefficient of determination between the predicted DMY from automated and manually derived polygons in the late-spring2018_2 was moderate (R² = 0.44). This could be due to smaller plant height architecture at that particular time, which coincided with early drought onset in the spring season of 2018. The canopy of these smaller plants may have reflected weak sound echoes from the surface of each plant, which may result in a greater variability in height values that resulted lower R².

Biomass yield distribution of perennial ryegrass cultivars may vary across growing seasons and years (Brock et al., 1996; Chapman et al., 2015). These variations mainly occur due to genetic merit of productivity and seasonal growth rate differences of cultivars during the vegetative growth stage. The ability to measure how much differences exist mainly depends on destructive harvest or visual score, which is restricted by low-throughput (Walter et al., 2012; Ghamkhar et al., 2019). Sensor-based non-destructive biomass yield predictive methods have the potential to replace this with applicability in large scale screening of breeding lines across seasons and years. Generally, individual plants NDVI, plant height and predicted biomass estimations showed variations across seasons. The distribution of NDVI values was wider in the winter2018 than late spring2018_1 and late spring2018_2. This indicates the various genotypes responded differently to cold season conditions (Förster et al., 2018). However, plant height and DMY predictions in the late spring2018_1 indicated higher values in which ryegrass appeared to have rapid growth and faster accumulation of dry matter in the spring season (Wang et al., 2019).

In this study, NDVI, plant height and predicted biomass estimations from the UAS and PhenoRover showed obvious temporal and spatial variabilities that mainly arises from genetic and environmental interactions. In the winter2018 harvest, NDVI histograms of some cultivars (Cultivar ID 8, 10, 13, 25, 45, 46, 47, and 48) indicated they were at NDVI saturation point before the harvest. This was further confirmed by the estimates of the plant height and DMY of these cultivars. These methods of accurate non-destructive estimates may provide a better perspective to facilitate breeding selection and ranking of ryegrass genotypes. The application of this non-destructive biomass prediction will enable the development of seasonal forage value indices (FVI) which combine biomass yield and economic values to select best performing commercial cultivars (Leddin et al., 2018; Giri et al., 2019). At the late-spring2018_1 harvest, plant height increased, and NDVI values saturated for all cultivars. The increase in plant height was from the fast growth rate in the spring season and NDVI value peaked and reached saturation before harvest. This result agrees with previous studies in wheat (Li et al., 2019) and buffelgrass (Mutanga and Skidmore, 2004) where NDVI tended to saturate at high canopy density. However, in the late-spring2018_2 lower NDVI, plant height and predicted biomass yields were observed for some plants. This indicates there was an early-onset of a late spring drought stress, that can restrict leaf growth (Cyriac et al., 2018), and plants with a thin canopy density are expected to reflect limited echoes and NDVI values (Condorelli et al., 2018; Yuan et al., 2018).

Experimental plots in this study were represented by three sets of rows containing 32 spaced-planted individuals each (total 96 plants/plot), in which individual plants level biomass prediction has reflected an accurate representation of plots biomass prediction. The results showed a significantly higher coefficients of determination between the predicted and measured DMY across three seasons with R² = 0.19, 0.75, and 0.81 at the winter2018, late-spring2018_1 and late-spring2018_2 harvests, respectively, and RMSE values ranging from 0.09 to 0.21 kg/plot. Subsequently, the overall regression of the three harvests indicated significantly high coefficients of determination (R² = 0.94) and RMSE 0.21 kg/plant. Similar coefficients of determination trends (R² = 0.76-0.89) were observed in previously reported pasture studies when various sensors were used to estimate biomass yield of tall fescue and ryegrass (Safari et al., 2016; Ghamkhar et al., 2019).

In this study, we validated a computational workflow for image acquisition, processing, and analysis to predict the biomass yield based on vegetative indices and plant height measurements of 48,000 ryegrass plants. Previous studies indicated the use of NDVI for ranking cultivars of ryegrass (Wang et al., 2019), field pea, canola, and spring wheat grain yield (Brian McConkey et al., 2004) and lint yield in cotton (Hugie et al., 2018). Considering the plant height and NDVI as a surrogate to predict DMY of individual and plot-level plants, there is a great potential to apply our workflow to be used for ranking of genotypes and cultivars across growing seasons and years. Our method also has the potential to further pave a foundation for the application of sensor technologies combined with GS that leads to faster rates of genetic gain in forages. To further improve the stability and accuracy of the automated workflow, other traits including vegetative indices other than NDVI, volume and leaf area features may be included by incorporating measurements through sensors including LiDAR and various camera systems (Ghamkhar et al., 2019; Yates et al., 2019).

Conclusion

The computational phenotyping workflow developed in this study indicated that automatically built polygons could replace manually derived polygons for plant height and NDVI accurate extraction from 48,000 plants. The NDVI, plant height, predicted and measured DMY showed accurate seasonal and spatial variations, which may contribute to identifying differences among genotypes of varying genetic background. Thus, the proposed computational phenotyping workflow proved to enable non-destructive and objective prediction of individual plants NDVI, plant height and DMY in three different seasons. Subsequently, the observed plot level linear relationship between predicted and measured DMY indicated capability of the developed computational workflow to estimate plot-level DMY from the individual plants’ prediction. Therefore, this study demonstrated the potential of the validated computational workflow to conduct seasonal and annual ranking of breeding lines and commercial cultivars that may be used for the application of genomic selection that paves to greater rates of genetic gain in forages.

References

• 1

  BarrachinaM.CristóbalJ.TullaA. F. (2015). Estimating above-ground biomass on mountain meadows and pastures through remote sensing.Intern. J. Appl. Earth Observ. Geoinform.38184–192. 10.1016/j.jag.2014.12.002

  • CrossRef
  • Google Scholar

• 2

  BarrettB. A.FavilleM. J.NicholsS. N.SimpsonW. R.BryanG. T.ConnerA. J. (2015). Breaking through the feed barrier: options for improving forage genetics.Anim. Product. Sci.55883–892.

  • Google Scholar

• 3

  BendigJ.BoltenA.BennertzS.BroscheitJ.EichfussS.BarethG. (2014). Estimating biomass of barley using crop surface models (CSMs) derived from UAV-based RGB imaging.Remote Sens.610395–10412. 10.3390/rs61110395

  • CrossRef
  • Google Scholar

• 4

  Borra-SerranoI.De SwaefT.MuylleH.NuyttensD.VangeyteJ.MertensK.et al (2019). Canopy height measurements and non-destructive biomass estimation of Lolium perenne swards using UAV imagery.Grass Forage Sci.74356–369.

  • Google Scholar

• 5

  Brian McConkeyP. B.LafondG. P.MoulinA.PelcatY. (2004). Optimal time for remote sensing to relate to crop grain yield on the Canadian prairies.Can. J. Plant Sci.8497–103. 10.4141/p03-070

  • CrossRef
  • Google Scholar

• 6

  BrockJ.HumeD.FletcherR. (1996). Seasonal variation in the morphology of perennial ryegrass (Lolium perenne) and cocksfoot (Dactylis glomerata) plants and populations in pastures under intensive sheep grazing.J. Agric. Sci.12637–51. 10.1017/s0021859600088791

  • CrossRef
  • Google Scholar

• 7

  BrocksS.BarethG. (2018). Estimating barley biomass with crop surface models from oblique RGB imagery.Remote Sens.10:268. 10.3390/rs10020268

  • CrossRef
  • Google Scholar

• 8

  BusemeyerL.RuckelshausenA.MöllerK.MelchingerA. E.AlheitK. V.MaurerH. P.et al (2013). Precision phenotyping of biomass accumulation in triticale reveals temporal genetic patterns of regulation.Sci. Rep.3:2442.

  • Google Scholar

• 9

  ChapmanD.MuirP.FavilleM. (2015). Persistence of dry matter yield among New Zealand perennial ryegrass (Lolium perenne L.) cultivars: insights from a long-term data set.J. New Zeal. Grassl.77177–184.

  • Google Scholar

• 10

  ChawadeA.van HamJ.BlomquistH.BaggeO.AlexanderssonE.OrtizR. (2019). High-throughput field-phenotyping tools for plant breeding and precision agriculture.Agronomy9:258. 10.3390/agronomy9050258

  • CrossRef
  • Google Scholar

• 11

  CondorelliG. E.MaccaferriM.NewcombM.Andrade-SanchezP.WhiteJ. W.FrenchA. N.et al (2018). Comparative aerial and ground based high throughput phenotyping for the genetic dissection of NDVI as a proxy for drought adaptive traits in durum wheat.Front. Plant Sci.9:893. 10.3389/fpls.2018.00893

  • CrossRef
  • Google Scholar

• 12

  CyriacD.HofmannR. W.StewartA.SathishP.WinefieldC. S.MootD. J. (2018). Intraspecific differences in long-term drought tolerance in perennial ryegrass.PLoS One13:e0194977. 10.1371/journal.pone.0194977

  • CrossRef
  • Google Scholar

• 13

  DeviaC. A.RojasJ. P.PetroE.MartinezC.MondragonI. F.PatinoD.et al (2019). High-throughput biomass estimation in rice crops using UAV multispectral imagery.J. Intellig. Robot. Syst.96573–589. 10.1007/s10846-019-01001-5

  • CrossRef
  • Google Scholar

• 14

  DobbelsA. A.LorenzA. J. (2019). Soybean iron deficiency chlorosis high throughput phenotyping using an unmanned aircraft system.Plant Methods15:97.

  • Google Scholar

• 15

  FörsterL.GrantJ.MichelT.NgC.BarthS. (2018). Growth under cold conditions in a wide perennial ryegrass panel is under tight physiological control.PeerJ6:e5520. 10.7717/peerj.5520

  • CrossRef
  • Google Scholar

• 16

  GebremedhinA.BadenhorstP.WangJ.GiriK.SpangenbergG.SmithK. (2019a). Development and validation of a model to combine NDVI and plant height for high-throughput phenotyping of herbage yield in a perennial ryegrass breeding program.Remote Sens.11:2494. 10.3390/rs11212494

  • CrossRef
  • Google Scholar

• 17

  GebremedhinA.BadenhorstP. E.WangJ.SpangenbergG. C.SmithK. F. (2019b). Prospects for measurement of dry matter yield in forage breeding programs using sensor technologies.Agronomy9:65. 10.3390/agronomy9020065

  • CrossRef
  • Google Scholar

• 18

  GhamkharK.IrieK.HagedornM.HsiaoJ.FourieJ.GebbieS.et al (2019). Real-time, non-destructive and in-field foliage yield and growth rate measurement in perennial ryegrass (Lolium perenne L.).Plant Methods15:72.

  • Google Scholar

• 19

  GiriK.ChiaK.ChandraS.SmithK. F.LeddinC. M.HoC. K. M.et al (2019). Modelling and prediction of dry matter yield of perennial ryegrass cultivars sown in multi-environment multi-harvest trials in south-eastern Australia.Field Crops Res.243:107614. 10.1016/j.fcr.2019.107614

  • CrossRef
  • Google Scholar

• 20

  HaghighattalabA.González PérezL.MondalS.SinghD.SchinstockD.RutkoskiJ.et al (2016). Application of unmanned aerial systems for high throughput phenotyping of large wheat breeding nurseries.Plant Methods12:35.

  • Google Scholar

• 21

  HayesB. J.CoganN. O. I.PembletonL. W.GoddardM. E.WangJ.SpangenbergG. C.et al (2013). Prospects for genomic selection in forage plant species.Plant Breed.132133–143. 10.1111/pbr.12037

  • CrossRef
  • Google Scholar

• 22

  HerridgeR. P.MacknightR. C.BrownfieldL. R. (2018). Prospects for F1 hybrid production in ryegrass.N. Zeal. J. Agric. Res.621–11. 10.1080/00288233.2018.1559867

  • CrossRef
  • Google Scholar

• 23

  HugieK. L.BauerP. J.StoneK. C.BarnesE. M.JonesD. C.CampbellB. T. (2018). Improving the precision of NDVI estimates in upland cotton field trials.Plant Phenom. J.11–9. 10.2135/tppj2017.09.0009

  • CrossRef
  • Google Scholar

• 24

  HuntC. L.JonesC. S.HickeyM.HatierJ. H. B. (2015). Estimation in the field of individual perennial ryegrass plant position and dry matter production using a custom-made high-throughput image analysis tool.Crop Sci.552910–2917. 10.2135/cropsci2015.02.0125

  • CrossRef
  • Google Scholar

• 25

  InostrozaL.AcuñaH.MunozP.VásquezC.IbáñezJ.TapiaG.et al (2016). Using aerial images and canopy spectral reflectance for high-throughput phenotyping of white clover.Crop Sci.562629–2637. 10.2135/cropsci2016.03.0156

  • CrossRef
  • Google Scholar

• 26

  Jimenez-BerniJ. A.DeeryD. M.Rozas-LarraondoP.CondonA. T. G.RebetzkeG. J.JamesR. A.et al (2018). High throughput determination of plant height, ground cover, and above-ground biomass in wheat with LiDAR.Front. Plant Sci.9:237. 10.3389/fpls.2018.00237

  • CrossRef
  • Google Scholar

• 27

  KhanZ.Rahimi-EichiV.HaefeleS.GarnettT.MiklavcicS. J. (2018). Estimation of vegetation indices for high-throughput phenotyping of wheat using aerial imaging.Plant Methods14:20.

  • Google Scholar

• 28

  LeddinC. M.JacobsJ. L.SmithK. F.GiriK.MalcolmB.HoC. K. M. (2018). Development of a system to rank perennial ryegrass cultivars according to their economic value to dairy farm businesses in south-eastern Australia.Anim. Product. Sci.581552–1558.

  • Google Scholar

• 29

  LiJ.Veeranampalayam-SivakumarA.-N.BhattaM.GarstN. D.StollH.Stephen BaenzigerP.et al (2019). Principal variable selection to explain grain yield variation in winter wheat from features extracted from UAV imagery.Plant Methods15:123.

  • Google Scholar

• 30

  LinZ.CoganN. O.PembletonL. W.SpangenbergG. C.ForsterJ. W.HayesB. J.et al (2016). Genetic gain and inbreeding from genomic selection in a simulated commercial breeding program for perennial ryegrass.Plant Genom.9:46.

  • Google Scholar

• 31

  LootensP.RuttinkT.RohdeA.CombesD.BarreP.Roldán-RuizI. (2016). High-throughput phenotyping of lateral expansion and regrowth of spaced Lolium perenne plants using on-field image analysis.Plant Methods12:32.

  • Google Scholar

• 32

  McDonaghJ.O’DonovanM.McEvoyM.GillilandT. J. (2016). Genetic gain in perennial ryegrass (Lolium perenne) varieties 1973 to 2013.Euphytica212187–199. 10.1007/s10681-016-1754-7

  • CrossRef
  • Google Scholar

• 33

  MutangaO.SkidmoreA. K. (2004). Narrow band vegetation indices overcome the saturation problem in biomass estimation.Intern. J. Remote Sens.253999–4014. 10.1080/01431160310001654923

  • CrossRef
  • Google Scholar

• 34

  PunalekarS. M.VerhoefA.QuaifeT. L.HumphriesD.BerminghamL.ReynoldsC. K. (2018). Application of Sentinel-2A data for pasture biomass monitoring using a physically based radiative transfer model.Remote Sens. Environ.218207–220. 10.1016/j.rse.2018.09.028

  • CrossRef
  • Google Scholar

• 35

  RahmanM. M.LambD. W.StanleyJ. N.TrotterM. G. (2014). Use of proximal sensors to evaluate at the sub-paddock scale a pasture growth-rate model based on light-use efficiency.Crop Pasture Sci.65400–409.

  • Google Scholar

• 36

  Rueda-AyalaP. V.PeñaM. J.HöglindM.Bengochea-GuevaraM. J.AndújarD. (2019). Comparing UAV-based technologies and RGB-D reconstruction methods for plant height and biomass monitoring on Grass Ley.Sensors19:535. 10.3390/s19030535

  • CrossRef
  • Google Scholar

• 37

  SafariH.FrickeT.ReddersenB.MöckelT.WachendorfM. (2016). Comparing mobile and static assessment of biomass in heterogeneous grassland with a multi-sensor system.J. Sens. Sens. Syst.5301–312. 10.5194/jsss-5-301-2016

  • CrossRef
  • Google Scholar

• 38

  ShiY.ThomassonJ. A.MurrayS. C.PughN. A.RooneyW. L.ShafianS.et al (2016). Unmanned aerial vehicles for high-throughput phenotyping and agronomic research.PLoS One11:e0159781. 10.1371/journal.pone.0159781

  • CrossRef
  • Google Scholar

• 39

  SmithK. F.TasneemM.KearneyG. A.ReedK. F. M.LeonforteA. (2001). Evaluation of herbage yield in a forage grass breeding program: comparison of visual rating versus measurement in single-row plots or swards.Anim. Product. Sci.411161–1166.

  • Google Scholar

• 40

  TattarisM.ReynoldsM. P.ChapmanS. C. (2016). A direct comparison of remote sensing approaches for high-throughput phenotyping in plant breeding.Front. Plant Sci.7:1131. 10.3389/fpls.2018.001131

  • CrossRef
  • Google Scholar

• 41

  ViljanenN.HonkavaaraE.NäsiR.HakalaT.NiemeläinenO.KaivosojaJ. (2018). A novel machine learning method for estimating biomass of grass swards using a photogrammetric canopy height model, images and vegetation indices captured by a drone.Agriculture8:70. 10.3390/agriculture8050070

  • CrossRef
  • Google Scholar

• 42

  WalterA.StuderB.KöllikerR. (2012). Advanced phenotyping offers opportunities for improved breeding of forage and turf species.Ann. Bot.1101271–1279. 10.1093/aob/mcs026

  • CrossRef
  • Google Scholar

• 43

  WangJ.BadenhorstP.PhelanA.PembletonL.ShiF.CoganN.et al (2019). Using sensors and unmanned aircraft systems for high-throughput phenotyping of biomass in perennial ryegrass breeding trials.Front. Plant Sci.10:1381. 10.3389/fpls.2019.01381

  • CrossRef
  • Google Scholar

• 44

  YatesS.JaškūnėK.LiebischF.NagelmüllerS.KirchgessnerN.KöllikerR.et al (2019). Phenotyping a dynamic trait: leaf growth of perennial ryegrass under water limiting conditions.Front. Plant Sci.10:344. 10.3389/fpls.2019.00344

  • CrossRef
  • Google Scholar

• 45

  YuanW.LiJ.BhattaM.ShiY.BaenzigerS. P.GeY. (2018). Wheat height estimation using LiDAR in comparison to ultrasonic sensor and UAS.Sensors18:3731. 10.3390/s18113731

  • CrossRef
  • Google Scholar

• 46

  YueJ.YangG.LiC.LiZ.WangY.FengH.et al (2017). Estimation of winter wheat above-ground biomass using unmanned aerial vehicle-based snapshot hyperspectral sensor and crop height improved models.Remote Sens.9:708. 10.3390/rs9070708

  • CrossRef
  • Google Scholar