High‐Throughput Approaches for Phenotyping Alfalfa Germplasm under Abiotic Stress in the Field

Core Ideas Remote sensing technologies enable rapid and nondestructive phenotyping of plants in the field. Biomass estimate accuracy from images and sensors was similar to yields harvested manually. Biomass yield variation identified the most productive accessions under low‐input conditions.

High-throughput phenotyping technologies enable monitoring plant growth and development nondestructively throughout the growing season. Crop losses associated with abiotic and biotic factors threaten the sustainability of crops used for feed, fiber, and fuel. The process to develop improved cultivars with enhanced crop yields includes phenotyping hundreds of plants under a target set of conditions. However, the manual collection of data is often laborious and time consuming. Strategies that integrate remote sensing technologies including unmanned aerial vehicles (UAVs) and sensors mounted on "phenomobiles" can streamline phenotyping efforts in plant breeding programs. The objectives of this study were to compare the phenotypic data collected from the field using UAVs, sensors, and manual approaches and to assess the potential of high-throughput approaches to rank the productivity of alfalfa (Medicago sativa L.) accessions growing in the field. Phenotypic data were collected from 100 alfalfa accessions established and grown under low-input conditions. Specific traits evaluated using both UAVs and sensors mounted on a phenomobile prior to physically harvesting the biomass during four harvests include biomass yield, plant height, normalized difference vegetation index, leaf area index, and ground coverage. The results from both the UAV and the sensors were highly correlated to the physical measurements obtained for the multiple traits evaluated. Therefore, field-based high-throughput phenotyping strategies represent a viable option for efficiently screening germplasm in the field to increase phenotyping efficiencies in plant breeding programs. I ncreased demand for agricultural products, erratic rainfall patterns, emerging pests and pathogens, and competing land-use interests are some of the driving forces behind the need to develop enhanced cultivars (Tilman et al., 2011). Additionally, competing interests between agriculture, urban development, and other uses will require the development of enhanced cultivars that are able to grow and thrive under less than optimal conditions, oftentimes in marginal soils, and be more productive with fewer inputs including fertilizers. Soils deficient in N and P, key elements for plant growth, are common in the United States (Hinsinger, 2001;White and Hammond, 2008;Smith et al., 2014) and other parts of the world. Nutrient deficiencies or toxicities are more apparent in acidic soils, causing discoloration of leaves due to the loss of chlorophyll, thus leading to decreased photosynthetic ability and lower yields. Fertilizers are often used to mitigate yield gaps, but this process increases production costs (Mikkelsen, 2004). Additionally, many crops are inefficient at P uptake and therefore only between 20 and 30% of P fertilizers are actively used by plants (López-Arredondo et al., 2014), thus contributing to fertilizer runoff with negative environmental impacts. Opportunities to identify natural variation for efficient biomass conversion under limited soil resources could result in cultivars with similar yields under reduced on-farm inputs.
Plant genetic resources including released cultivars, landraces, wild or related species, elite germplasm or breeding lines, and genetic stocks are evaluated to identify variation for traits of interest (yield, stress tolerance, or disease resistance) captured within the beneficial alleles of certain genes (Haussmann et al., 2004). Germplasm collections are available for a diverse range of agricultural species and can be used to identify and expand the genetic variation for improving agricultural crops to alleviate some of the genetic bottlenecks reported in certain cultivated species (Tanksley and McCouch, 1997;van de Wouw et al., 2009).
Alfalfa is a perennial legume species capable of symbiotic N 2 fixation. This legume is widely grown in North America, South America, and Europe, and these countries combined account for approximately 90% of global production (Michaud et al., 1988). The aboveground biomass is a high-quality forage consumed directly through grazing or harvested as hay multiple times throughout the growing season. Reproducible and accurate plant phenotyping is a critical component to develop enhanced cultivars based on the identification of genetic variation for key traits.
Thousands of accessions from the alfalfa germplasm collection were characterized to identify a subset of 200 core accessions representing the diversity of the collection based on their phenotypic characteristics (Basigalup et al., 1995). However, the phenotyping process for diverse plants in the field can be labor intensive and time consuming (Montes et al., 2011;Fahlgren et al., 2015) and only documents a data point at a specific moment. Strategies that would enable data collection nondestructively throughout the growing season would generate more robust data, for example to enable breeders to distinguish plants with early vigor and lower yields later in the season in contrast to plants with slower establishment but higher biomass yield at the end of the season. However, the sole collection of nondestructive data considerably limits the phenotyping of plants in field. Indeed, visual rating of plant growth or disease impacts is well known to be very subjective, depending on the experimenter, and would also be time consuming. The survivability of individuals under low-input conditions would be determined by differences existing between accessions. Further, capturing the variability between plants that merely survive and the surviving plants that also have high yields under stress are key distinctions to identify the best parents to use for breeding.
The emerging availability of remote sensing technologies including high-throughput phenotyping platforms can facilitate screening of thousands of plants under different conditions (White et al., 2012;Cowley et al., 2014;Rembold et al., 2015). Systems that use sensors for data acquisition under varying payloads are available (Kirchgessner et al., 2017;Virlet et al., 2017). However, these systems require multiple hours to perform a full cycle of data collection and, due to their fixed location, are limited to data obtained only within a certain plot area. Further, ground-based "phenomobiles" (Deery et al., 2014;Pittman et al., 2015Pittman et al., , 2016 and UAVs or "polycopters" (Hogan et al., 2017) are mobile and more versatile in terms of the number of plots, the types of sensors used and their ability to overcome uneven terrain or fences in grazing trials. Remote sensing technologies including UAVs and sensors can gather field data at high spatial resolution and temporal frequency. However, these technologies are subject to a lower payload capacity (maximum 2 kg for some UAVs) (Sankaran et al., 2015;Hogan et al., 2017), limited battery life, or fuel capacity and adverse weather conditions (wind, rain, or temperature extremes) may hinder their ability to operate.
Traits previously evaluated using sensor technologies include laser and sonic height, canopy temperature, and the normalized difference vegetation index (NDVI), among others. The NDVI is calculated using wavelengths within the near-infrared and visible regions of the electromagnetic spectrum (Tucker, 1979). The NDVI relates to chlorophyll content based on its absorption properties that can reflect the photosynthetic machinery of plants (Tattaris et al., 2016). Different types of cameras can be used to capture images based on the red-green-blue (RGB) spectrum, NDVI, or thermal images (Pittman et al., 2015(Pittman et al., , 2016Hogan et al., 2017;Kirchgessner et al., 2017;Virlet et al., 2017).
The successful deployment of high-throughput phenotyping strategies in breeding programs aimed at developing enhanced cultivars with efficient biomass conversion under low-input conditions should be developed and evaluated against conventional methods for phenotypic data collection. Therefore, the objectives of this research were to deploy remote sensing and manual data collection approaches to discriminate between alfalfa accessions for field-based performance, productivity, and biomass yield.

Field Site
On 28 Sept. 2015, 14 soil cores (4.44-cm diameter and 122 cm deep) were collected as duplicates from a field site located at the Noble Research Institute's Research Park Farm in Ardmore, OK, using a Giddings soil-sampling machine (Model GSRTS, serial no. 1397) mounted to a John Deere 6615 tractor. The soil cores were collected using a zero-contamination probe fitted with a plastic sleeve. Each soil core sample was subdivided into five sections corresponding to the different soil depths of 0 to 15, 15.1 to 30, 30.1 to 60, 60.1 to 91.5, and 91.5 to 122 cm. Soil cores were submitted for pH and nutrient composition (N, K, and P) analysis to ServiTech (Amarillo, TX). The information from the soil analysis enabled assignment of four of the 12 sampled parcels (Supplemental Table  S1) within the total field area to the four replications of a randomized complete block design to evaluate the alfalfa accessions. The parcels used in this study were selected based on the uniformity of the soil pH and consistent low P in four adjacent plot regions. The GPS coordinates of the field plot use for planting are: 34°11¢39.984² to 34°11¢38.76² N, 97°5¢1.32² to 97°4¢59.52² W, and an elevation of 262.6 m above sea level.

Plant Materials
Seeds from a total of 94 M. sativa plant introductions (PIs) and the check cultivar Vernal were obtained from the Germplasm Resources Information Network (GRIN, https://npgsweb.arsgrin.gov/gringlobal/search.aspx). The accessions evaluated were identified based on descriptors suggesting nutrient use efficiency, drought tolerance, or source of origin in regions known to have acidic and/or nutrient-deficient soils originating from multiple countries (Table 1 and Supplemental Table S2). PI 525455 is recommended as a genetic source to increase P concentration in alfalfa cultivars (Melton et al., 1989). Five entries corresponding to the named cultivars Bulldog 805 (Bouton et al., 1997) and Mariner IV, FSG420LH, FSG423ST, and FSG639ST (Allied Seed, contrasting for root growth) completed the list of 100 accessions. For each accession, 100 seeds were manually scarified using 180 grit sandpaper and planted in 11.43-cm pots containing Metromix professional soil (Sun Gro Horticulture) previously saturated to capacity with water on 10 Aug. 2016. The Metromix 360 soil contains small amounts of slow-release N. Seedlings were grown in the greenhouse for 5 d, and the 72 most uniform seedlings from each accession were transplanted into individual 96-cell flats with Metromix 360 soil also previously saturated with water; the seedlings were then grown in the greenhouse. Twelve days after germination, the plants were inoculated with a mix of Rhizobium (3,100,000 colony-forming units [cfu] g −1 for multiple plant species), Azotobacter (2,400,000 cfu g −1 ), and Bacillus (1,000,000 cfu g −1 for multiple plant species) (Micronoc Seeds Inoculant, Sono AG Products). Seedling plant growth was cut back to the third node on 3 Oct. 2016 to minimize variations associated with differences in seedling emergence and vigor. The experimental design was a randomized complete block with four replications and 10 individuals per replication for a total of 40 plants evaluated per accession as previously described for alfalfa spaced plants (Bhattarai et al., 2018). All plants were transplanted to the field site on 18 Oct. 2016. Each row plot contained 10 plants planted at 30.48-cm spacing between individual plants and 91.44-cm spacing between rows. Additionally, a border row of the cultivar Bulldog 805 was planted around the periphery of the entire plot area to minimize border effects. Low-input management practices were implemented, meaning that no fertilizers, herbicides, or insecticides were applied. Additionally, the field was rainfed (no irrigation), with a total 87.3 cm of precipitation for the duration of the trial (data from the Oklahoma Mesonet weather station).

Manual Data Collection, Storage, and Processing
After establishment, the number of surviving plants per accession were determined based on visual ratings considering the overall health and status of the plant. Relative growth was determined on 30 Nov. 2016 using a visual scale ranging from 1 (lowest) to 10 (highest) considering the average size and biomass of the 10 plants per row plot as a way to nondestructively detect differences among accessions by walking throughout the field plot area. The Table 1. One hundred alfalfa (Medicago sativa L.) accessions evaluated in the field. The Entry ID ranges from 1 to 100, and the corresponding Plant Introduction (PI) number was obtained from the Germplasm Resource Information Network (GRIN). The entries 1 and 96 to 100 correspond to commercial alfalfa cultivars that were used as checks. average visual score of the 10 individual plants does not capture the homogeneity or variation in plant growth among the 10 plants of the plot. The raw data were collected directly in the field using a tablet (Lenovo ThinkPad) containing the field plot map. The aboveground biomass was harvested manually at 5-cm stubble height four times-on 15 May, 13 June, 14 July, and 17 Aug, 2017-to determine fresh biomass (kg) per plot onsite as previously described (Bhattarai et al., 2018), using an Ohaus Defender 3000 scale (Parsippany, NJ) connected directly to a rugged tablet (Panasonic Toughpad). Samples placed inside labeled paper bags were dried inside forced-air ovens at 63°C for 2 d to determine shoot dry biomass. The number of plants remaining in each plot during each harvest was used to determine the fresh and dry shoot biomass per plant within the plot.

UAV Data Collection and Processing
A DJI Inspire 1 UAV equipped with a DJI Zenmuse X5 visual band sensor camera (resolution 0.6-1.1 cm per pixel for RGB images) and a Sentera single sensor multispectral camera (resolution 2.5 cm per pixel for NDVI images) (Supplemental Fig. S1A) was deployed for data collection. The total payload of this UAV was 560 g with both cameras, with an approximate 20-min flight time on a single battery charge. Images from the UAV were collected at an altitude of 29 m flown at solar noon to minimize the impacts from shadows during the image acquisition process and flown only when the wind speed was below 32 km h −1 . An operator with a remote pilot certification from the Federal Aviation Administration (FAA) flew the UAV in compliance with all regulations. The individual UAV-generated images were stitched together using AgiSoft Photoscan Professional Version 1.2.2 software (Agisoft LLC) to create orthomosaic images corresponding to the entire field site. The stitching process involves aligning photos and building a dense point cloud, mesh surface, and orthomosaic. During this process, individual photos were aligned and stitched with high accuracy and high image quality parameters. Finally, the orthomosaic images were geo-referenced using ArcGIS ArcMap software Version 10.3.1 (ESRI) by integrating the GPS coordinates from five ground control points positioned both at the center and in each corner of the field site (Supplemental Fig. S1). The corresponding coordinates for each of the five ground control points were determined using a Geo 7X handheld data collector (Trimble Geo Explorer).
Two models were developed using the ArcGIS model builder software to calculate and extract data from the multispectral and RGB orthomosaic images to estimate plant ground coverage (Fig. 1). These models were automated and comprised multiple ArcGIS tools to extract analytical information from the orthomosaic images. The NDVI model was used to extract the NDVI data from the UAV-generated multispectral images and can be used to represent the plant growth and vigor by comparing the spectral reflectance in the red (wavelength 575-700 nm) and near infrared (wavelength 830-870 nm) regions of the spectrum. The NDVI is the most common vegetation index to measure plant health based on the accumulation of chlorophyll in the leaves, was first reported by Rouse et al. (1974), and was calculated based on the standard NDVI equation and populated using the NDVI orthomosaic images: Delimiting the field plot area from these images involved the use of a mask outlining the perimeter of the field plot to remove roads, people, vehicles, or any other elements that could bias the data. The NDVI formula was applied to the image to create a heat map representing the numerical NDVI values ranging from 0 to 1, where values closer to 0 represent less healthy plants or bare ground and values closer to 1 represent healthier plants. A filtering step was applied to the heat map to remove the background effects or noise in the orthomosaic images prior to determining the plant NDVI values. To achieve this step, bareground NDVI values and plant NDVI values were segmented out using ArcGIS software tools. Any NDVI values £0.2 were considered bare-ground values, and any values higher than the bare-ground values were considered part of the NDVI values. The plant ground-coverage model was used to estimate the area of the ground covered by plants. This utilizes the multispectral images generated as described above. The field plot area was also extracted from the RGB orthomosaic using the image masking function in ArcGIS. The bare-ground background was then suppressed in the image, leaving the remaining pixels in the image corresponding to the values associated with plants only. Initially, a fishnet grid including the plant number and replicate number was created using the ArcGIS software to match the data from each row and column that was extracted from the collected image. Using the Zonal Statistic tool in ArcGIS, the fishnet grid was then applied to the masked image to extract the number of pixels as well as the total area covered by these pixels for each plot row directly into a spreadsheet. The total pixel area corresponds to the plant ground coverage and was used for further analysis.

Sensors (Phenomobiles) Data Collection and Processing
The spider refers to a gasoline-powered high-clearance tractor (Lee Spider, LeeAgra) mounted with an array of sensors (Supplemental Fig. S2) and used for field-based data collection as previously described (Pittman et al., 2016). The spider's hydraulic, front-end platform contains two sensing platforms with identical components. These platforms are fitted with two time-of-flight light lasers (Micro-epsilon Model ILR 1030) and one 120-MHz ultrasonic distance sensor (Senix TSPC 30S1-232) to determine plant height. Three active multispectral canopy sensors (Crop Circle models ACS-430, ACS-470, and DAS43X, Holland Scientific) documented ambient and canopy temperature, relative humidity, atmospheric pressure, photosynthetic active radiation, vegetation indices (normalized difference vegetation index [NDVI], inverted ratio vegetation index [IRVI], normalized difference red edge [NDRE], leaf area index vegetation index [LAI.VI], and canopy chlorophyll content vegetation index [CCC.VI]), and reflectance information. One Green Seeker (Trimble) equipped each sensor platform used to determine the NDVI values using the same NDVI formula shown in the previous section. This sensor emits bursts of red and infrared light and then measures the amount of each light wave reflected back from the plant. A global navigation satellite system (Trimble Omnistar NovaTel GPS Flexpak 6, 5-cm accuracy/error) enabled mapping of the values collected from the sensors to a particular position in the field and hence a specific row plot within the field area. Data were obtained at a rate of 10 Hz (for example, 10 samples s −1 ) and connected via USB to a field tablet (Panasonic Toughpad FZ-G1) and collected in real time utilizing the AgriLogger software (custom software application developed at the Noble Research Institute). The AgriLogger software served to retrieve, fuse, and write data from each instrument in real time to a .csv file format. The data were processed and extracted on a per-plot basis utilizing tools and applications in ArcGIS ArcMap 10.3 software (ESRI). Final data deliverables included per-plot summary statistics for all plant parameters. In addition, empirical Bayesian kriging surface models in ArcGIS ArcMap 10.3 were produced to generate heat maps for each of the plant growth parameters (Fig. 2B-G) similar to the ones generated from the UAV images. These surface models used the data points provided by the spider while driving through the field. Previously described models developed to estimate plant biomass (Pittman et al., 2015) and crude protein content (Pittman et al., 2016) in alfalfa were implemented in this study as described above. The frequency of data collection established was 2-wk intervals between 7 Dec. 2016 and 14 Sept. 2017, occasionally adjusted for weather conditions.

Calculation of Phenotypic Data Descriptors
Biomass is a parameter estimated using NDVI and plant height data from the spider only, as previously described (Pittman et al., 2015). The NDVI and ground coverage datasets were extracted following the use of the models described above that were highly correlated to the physical biomass measured in the field. The aboveground dry biomass was therefore estimated using the NDVI and the ground coverage dry weight (DW) from the UAV images as The statistical regression between the estimated dry biomass and the dry biomass collected during harvests was R 2 = 0.868 (in this study). This estimation was performed using the partial least square regression method, which can fit a regression model when the correlation between the predicator variables is significant (multi-collinearity problems). The model considered the physical dry weight measurements from all four harvests, UAV-obtained NDVI, and ground coverage data. The normality assumption was verified by checking the distribution of Studentized residuals as well as checking the departure from the theoretical line using a Q-Q plot (data not shown).
The leaf area index (LAI), noted Spider_LAI.VI in this study, was collected using the sensors from the spider. To compare with variables coming from the UAV, several other formulas previously described were used to estimate the LAI based on plant height, coverage, or NDVI. For example, the equation of Thorp et al. (2015) integrates plant height and coverage and was used to analyze the data collected in this study. This equation combines the plant height data from the spider and the ground coverage data from the UAV: where b = 5.5, being a constant described by Thorp et al. (2015), c and h are the canopy ground coverage and the plant height, respectively, and c max and h max are the maximum ground coverage and plant height expected during the growing season. Finally, the plant volume estimate was based on the plant height values obtained from the spider and the plant ground coverage obtained from the UAV.
An analysis of variance (ANOVA) explored the sources of variation due to genotype, replication, and harvest and the interactions between these including genotype ´ replicate, evaluated using ProcGLM (a general linear model procedure) in the SAS 9.4 software (SAS Institute). A split-plot design with a randomized complete block design was used to analyze the relationships between factors and response variables. The four harvests were used as a block effect. Replicates (n = 4) were assigned to the main plots with each block (harvest). The genotypes (entries) were used as subplot levels and then assigned within each main plot. The interaction between harvests and replicates was used as the error term of the main plot instead of the residual mean squared error, which was, therefore, the error term from the subplot. The ProcGLM with two test statements was implemented to generate the ANOVA table for each response variable.
Analysis of variance was used to determine effects due to genotype for each of the traits evaluated. The k-means clustering method was used as an alternative way of replacing the previous multiple comparison. The k-means clustering partitions observations into k clusters based on the "distance" between observations, which is similar to multiple comparison. The optimal number of k clusters was determined separately for each harvest and for all four harvests combined using the pamk function in the R package fpc.
Pearson correlation coefficients within each harvest obtained from manual harvesting of the plants including fresh and dry weight and the multiple parameters that were collected using the UAV and spider were determined using custom R scripts (R Development Core Team, 2008). Correlation between variables was determined based on p values <0.001, and correlation values >99.2 were rounded up to 100 in the correlation heat maps. The practical value of the correlation heat maps is to identify phenotypic trait(s) collected using remote sensing (RS) technologies that could be used to predict the biomass collected manually.

Field Site Conditions
The soil pH in the field area ranged from pH 5.1 to 5.5 in the first 30 cm of soil, and the P content of the upper layer of the soil ranged from 9 to 15.7 kg ha −1 (Supplemental Table S1). A more neutral pH of 6 to 7 and a P range between 67 and 73 kg ha −1 are considered optimal for alfalfa growth (Kelling, 2000;Undersander et al., 2011;Oklahoma Cooperative Extension Service, 2017), thus confirming the less than optimal pH and P-deficient conditions in the field. The two other nutrients to establish alfalfa and contribute to the longevity of the stand include N and K, and these should range from 22 to 33 and 175 to 437 kg ha −1 , respectively (Kelling, 2000;Undersander et al., 2011;Oklahoma Cooperative Extension Service, 2017). The N levels were within this range in the first 15 cm of the soil but decreased drastically in the deeper layers of soil. Inoculated alfalfa seedlings should be capable of symbiotic N 2 fixation. The levels of K in the soil were also lower than the recommended amounts and provided a test environment that included nutrient-deficient soil in which to evaluate alfalfa productivity with no inputs.

Persistence of Alfalfa Accessions during the Growing Season
Genetic variation for local adaptability and winter survival were identified in the alfalfa plants evaluated. The total number of surviving plants counted manually declined from 4000 to approximately 2700 plants in 12 mo (planting in October 2016 to October 2017 when fall dormancy ratings were obtained) (Fig. 3). Despite the strong selection pressure from the nutrient-deficient soils, 42 accessions had >31 plants across the four replications survive during the entire growing season, and 16 of them had 36 or greater surviving plants (Fig. 3). For example, all 40 plants for Entry 23 (PI 511303) persisted throughout the entire growing season. In contrast, some accessions did lose >75% of the plants-probably the result of the absence of cold tolerance in these accessions given the number of surviving plants surveyed in the spring. For example, two accessions had fewer than five survivors and three accessions had fewer than 10 survivors. The remaining accessions, with about 50 and 25% plant loss overall, also saw most of the decline in surviving plants during the spring, with another decline in survival after the fourth harvest (Fig. 3). All accessions, regardless of the number of surviving plants, were included in the data analysis.

Biomass Yield and Other Agronomic Traits
Alfalfa accessions had a range in total aboveground dry weight biomass (collected destructively) when evaluated under low-input field conditions between 0.04 to 1.3 kg per plot (Table 1; Fig. 4A). The range in dry biomass for each harvest was between almost 0 to 0.560 kg for Harvest 1, 0.016 to 0.384 kg for Harvest 2, 0.012 to 0.246 kg for Harvest 3, and 0.001 to 0.228 kg for Harvest 4 (Fig.  4A). A harvest effect was identified at p < 0.001 (data not shown), with the highest dry biomass obtained during Harvest 1 and a reduction in biomass on a per-harvest basis as the season progressed until the final Harvest 4 (Fig. 4A). The dry biomass on a per-harvest basis was highly correlated with the total dry biomass from the four harvests combined based on R 2 values ranging between 0.73 and 0.94 (Fig. 4B). Harvest 1 also had a wider spread in biomass production between the highest and lowest yielding accessions compared with a smaller range in Harvests 3 and 4 (Fig. 4C).
A genotype effect due to differences among accessions for biomass (fresh and dry weight), NDVI, and other traits was significant based on the ANOVA (Supplemental Table S3). The k-means clustering method ranked and grouped accessions based on their biomass yield per harvest ( Fig. 5; Supplemental Table S4) and combined for the total biomass production for all four harvests (Fig.  6). The number of clusters generated was determined separately for each harvest as well as the total combined yield throughout the growing season and varied across the four harvests (Fig. 5). Entries 57, 66, and 72, corresponding to PI 643396, PI 304220, and PI 434600, respectively, were consistently in the cluster with the highest biomass for all four harvests (Supplemental Table S4), significantly different from the others. Two of the three top accessions (PI 304220 and PI 434600) originated from Argentina. There was no correlation between the biomass (fresh and dry weight) and the number of surviving individuals. Indeed, some accessions with a low biomass had a higher number of surviving individuals than some accessions with fewer plants but higher biomass per plant.

Correlation Analysis between Manual Data and Remote Sensing Data
Heat map correlations between the data collected manually and data obtained using RS technologies (spider and UAV) reveal significant correlations for all four harvests (Fig. 7). The data collected using sensors mounted on the spider include NDVI, plant height (from laser and ultrasonic sensors), and canopy temperature. The first three variables were used within a model to estimate biomass and crude protein content. The data obtained from the UAV include NDVI (extracted from the multispectral images) and plant ground coverage (extracted from the RGB images) (Fig. 1). The shade of color within the correlation matrix shown in blue or red in Fig. 7 corresponds to the correlation between the two variables and are positive or negative, respectively. The numbers on the opposite diagonal reflect the values corresponding to the correlation coefficient r. Given the harvest effect, correlations for each harvest are shown separately. For all harvests (1-4), the freshweight and dry-weight biomass values were positively or negatively correlated (r values greater than ±0.50) with the sensor and UAV generated data (Fig. 7). The correlations for the NDVI ratings ranged between 0.66 to 0.81 from the spider data and from 0.48 to 0.68 for the UAV data. Differences in how the time-of-flight laser (single beam) and ultrasonic sensor (broad area) collect plant height data resulted in some differences in the plant height generated from the two systems mounted on the spider. The plant height data generated from the ultrasonic sensor and the laser were both highly correlated to the biomass data (r value range of 0.68-0.80 for the ultrasonic sensor and 0.65-0.82 for the laser).
The total biomass for each accession (estimated in tonnes ha −1 ) and the crude protein content (as a percentage of the biomass) were determined based on the NDVI data and both the laser and ultrasonic height data collected using the spider, as described by Pittman et al. (2015Pittman et al. ( , 2016. The estimated biomass based on the sensor data was positively correlated with the manually harvested dry biomass based on r values ranging between 0.71 and 0.85. The correlations between the measured biomass and the estimated crude protein ranged from 0.70 to 0.82 (Fig. 7).
The NDVI values and plant ground coverage data collected using the UAV were used to develop a model to estimate the biomass dry weight. The model-based values for dry weight were highly correlated to the manually harvested dry weight from all four harvests, and the correlations ranged between 0.88 and 0.91. Further, the correlation between the estimated biomass obtained from the sensors with the UAV ranged between 0.74 and 0.83 (Fig.  7), even though they were fully independent datasets. Differences Fig. 5. Clustering of the aboveground dry biomass for the accessions evaluated ranked by increasing biomass. Each color represents a different cluster based on the statistical k-means clustering. The clustering for each harvest was determined based on statistically significant differences (p < 0.001) for data within each harvest and thus can differ between harvests. Each number above the dot represents the entry number for each accession.
between the plant growth habit (such as between prostrate vs. erect growth) could result in different values for plant height even though they may have similar biomass production. Therefore, the model developed for plant coverage considers these differences in the image analysis. The biomass estimates generated from this model were the most highly correlated with the manually harvested dry biomass data based on correlation values between 0.88 and 0.92.
Data collected independently from both the sensors and UAV images were used to generate values for the leaf area index (LAI).
The LAI values obtained directly as output from the Crop Circle sensors mounted on the spider had correlation values between 0.60 and 0.78 with the manual harvest data. The UAV-derived NDVI values combined with the Saito et al. (2001) formula generated LAI values with correlation values between 0.58 and 0.81 for all harvests (Fig. 7). Combining the plant height data collected from the laser and ultrasonic sensors mounted on the spider with the plant ground cover obtained from the UAV images enabled estimation of both the volume of the plants and the LAI using the Thorp et al. (2015) formula. The correlation values between plant volume and the manually collected dry biomass ranged between 0.87 and 0.90 and between 0.85 and 0.89 for the LAI (Fig. 7). The correlation values among the three estimates of LAI generated independently ranged between 0.51 and 0.83 across all four harvests (Fig. 7).
The canopy chlorophyll content vegetation index, measured with sensors mounted on the spider, had a positive correlation with both fresh-weight and dry-weight biomass estimates, but the correlations were lower than for the ground coverage and the biomass estimates from the models developed (Fig. 7). Negative correlations between the canopy temperature and the collected biomass obtained during all four harvests ranged between −0.16 and −0.64 (Fig. 7). Therefore, plants with higher canopy temperatures during the summer months had lower biomass yields and vice versa.

Discussion
Breeders aim to develop improved cultivars that can tolerate a variety of abiotic stress conditions, and these efforts oftentimes include evaluating specific traits that allow plants to adapt to and be grown under varying nutrient and water deficiencies in the soil (Paez-Garcia et al., 2015). Opportunities to implement RS approaches can streamline alfalfa phenotyping efforts in the field, as they are less labor intensive and time consuming than manual phenotyping. Differences in the biomass yield and other traits evaluated are representative of the genetic variability among the alfalfa accessions evaluated. Alfalfa grows best on soils that are well drained, neutral in pH, and have high fertility (Sheaffer and Evers, 2007), and thus the low-input growing conditions in the field represented a multifactor stress environment for screening. Additional factors affecting the adaptability and persistence of many forage species such as alfalfa include the absence of cold tolerance, drought susceptibility, and diseases, among others (Sleper and Poehlman, 2006, p. 335-360). A reduction in the number of surviving alfalfa plants after the winter probably reflects their sensitivity to cold temperatures. Later in the growing season, disease susceptibility and stress from lack of irrigation were possible factors driving the reductions in productivity and stand losses. Individual plants with the capacity to tolerate low pH, deficient nutrient conditions, winter temperatures, diseases, and other stress factors would persist and continue to grow, while individuals sensitive to any of these factors would not survive. Further, the RS technologies used enabled the distinction between plants that merely survived vs. those that not only survived but also had higher yields. For example, some had a survivability of 75% but the plants were short and low yielding. Therefore, considering the number of survivors combined with the biomass represents a more viable approach for selecting individuals suitable to use for parents in breeding programs.
Accessions with higher yields than the check cultivars under these low-input conditions, as shown with the k-means method, were identified and thus represent a potential source of novel alleles to integrate into alfalfa breeding programs. Although plants have developed adaptive strategies to cope with nutrient limitations (such as modifying gene expression, altering carbon metabolism, exudation of organic acids, and enhancing nutrient use efficiency by modifying root structures) (Vance et al., 2003), the specific strategies deployed by high-yielding accessions that also enabled them to persist under these low-input conditions require further evaluation.
The development and use of high-throughput phenotyping technologies such as drones and sensors mounted on ground-based vehicles can be useful in multiple screening and breeding applications. These include evaluating diverse accessions from germplasm collections for allele mining and quantification of smaller differences in productivity among elite breeding lines (Furbank and Fig. 6. Clustering of the total aboveground dry biomass (in kg per plot) for the accessions evaluated ranked by increasing biomass from all four harvests and assigned to clusters based on statistically significant differences (p < 0.001).
Tester, 2011; Tattaris et al., 2016), monitoring insect pressure in crop species (Lamp, 2018), or predicting biomass yield. Zeng and Chen (2018) used near-infrared vegetation canopy reflectance measured in the field to estimate forage biomass yield and nutrient contents of a legume-grass mixture with an accuracy of prediction between 0.67 and 0.80. The sensor-based approaches described in this study had higher correlations with manually harvested data points. Previous studies have used the canopy temperature Fig. 7. Heat map correlations between the manual data collection (FW = fresh weight, DW = dry weight) and data obtained from the spider tractor and the unmanned aerial vehicle (UAV) for Harvests 1 to 4. NDVI, normalized difference vegetation index; IRVI, inverted ratio vegetation index; NDRE, normalized difference red edge; LAI.VI, leaf area index vegetation index; CCC.VI, canopy chlorophyll content vegetation index; RE, red edge color reflectance; NIR, near-infrared color reflectance; RED, red color reflectance. The numbers within the diagonal indicate the value of the positive (blue) or negative correlation (red). The software rounds up the correlation values that are ³99.2 to 100. or canopy reflectance to estimate alfalfa biomass yield (Noland et al., 2018). Comparisons of plant height estimates obtained from laser and ultrasonic sensors and physical measurements on alfalfa, bermudagrass (Cynodon dactylon L.) and wheat (Triticum aestivum L.) in the field indicated that the differences between methods was relatively small (about an average error of 11%) between predicted and harvested values (Pittman et al., 2015). Although individual metrics or a combination of plant health indices (NDVI, LAI, plant height) were highly correlated with the phenology estimates in alfalfa, the accuracy of sensor-based approaches for these and other traits can be affected by the canopy architecture and leaf shape unique to each species.
Factors such as row spacing, crop species, soil coverage, plant architecture, and ground coverage can affect the estimation of NDVI and LAI (Redowan and Kanan, 2012). Therefore, combining multiple sensor or image-based data points could increase the accuracy of the models used to estimate biomass yields. In this study, we evaluated the correlation of both individual data points and a combination of RS-based data to generate models for estimating biomass yields in alfalfa. Estimates of total biomass and crude protein were determined using models highly correlated with dry-weight biomass previously described (Pittman et al., 2015(Pittman et al., , 2016. Additionally, the dry-weight biomass values were also estimated by developing a new model using the NDVI data from the "false color" UAV images and the plant ground-coverage data obtained from the RGB images to replace the plant height used in the study by Pittman et al. (2016). The correlation values for dry-weight biomass using this combined approach increased to 0.90. Similarly, Thorp et al. (2015) combined the plant height from sensors mounted on phenomobiles with the plant ground-coverage data obtained from the UAV.
Remote sensing phenotyping approaches do have certain limitations. Some factors to consider include the difference in reflectance among crops grown as mixtures, the effect of the position of the sun, which could result in shadows in the field from plant canopies, and the light intensity during sunny vs. cloudy days. Many of these factors, however, are not unique to sensor-based approaches and these could be more significant during manual data collection during which variation in the sun's position between collection of the first data point and the last data point would be much broader. Further, pre-calibration, setting up GPS ground control points to geo-correct an image, running the spider or flying the UAVs, and developing prediction models are part of the process for data collection. Another consideration is the data processing, filtering, and stitching of UAV-collected images to assign a quantifiable value for a phenotypic trait to an individual plot or plant. As these processes become routine, the time required for each step decreases. Financial considerations can also play a role in the deployment of these RS technologies. Costs for the multiple sensors, cameras, phenomobiles, UAVs, FAA licenses, and related software continue to decrease with technological advances. Therefore, the development costs could be spread across multiple years and also shared by multiple breeding programs. Despite these factors, the correlations for RS phenotypic data and manually harvested biomass highlight the value of these approaches to reduce the time and labor needed to collect data.
Overall, multiple strategies for plant phenotyping complement each other with regard to manual harvesting. However, the availability of multiple UAV models at affordable price points would enable its rapid deployment by multiple programs for data collection in the field. While the accuracy may be slightly lower than for the spider, a suitable filtering data extraction model can mostly compensate. Further, deployment of UAVs requires fewer technical components and facilitates faster scanning of an area without the limitations of driving through the field, including partially flooded areas following a heavy rainfall, making determination of the height of the plants growing in the field feasible, regardless of how tall the plants are, and the method is not subject to the wear and tear on tires associated with using sensors mounted on tractors.
Once the infrastructure for RS approaches is in place, their successful deployment for phenotyping forage species for biomass yield and persistence as part of trait-development and breeding programs can be used to reliably identify differences among genotypes. Additional benefits of these and other sensor-based methods include relatively fast data acquisition and data processing, among others, as described by Deery et al. (2016). Further benefits will be obtained as these strategies are deployed to evaluate progenies resulting from multiple cycles of selection, intermating, and evaluation aimed at increasing yield and persistence under low-input conditions.