UAV-Based LiDAR for High-Throughput Determination of Plant Height and Above-Ground Biomass of the Bioenergy Grass Arundo donax

: Replacing fossil fuels with cellulosic biofuels is a valuable component of reducing the drivers of climate change. This leads to a requirement to develop more productive bioenergy crops, such as Arundo donax with the aim of increasing above-ground biomass (AGB). However, direct measurement of AGB is time consuming, destructive, and labor-intensive. Phenotyping of plant height and biomass production is a bottleneck in genomics- and phenomics-assisted breeding. Here, an unmanned aerial vehicle (UAV) for remote sensing equipped with light detection and ranging (LiDAR) was tested for remote plant height and biomass determination in A. donax. Experiments were conducted on three A. donax ecotypes grown in well-watered and moderate drought stress conditions. A novel UAV-LiDAR data collection and processing workﬂow produced a dense three-dimensional (3D) point cloud for crop height estimation through a normalized digital surface model (DSM) that acts as a crop height model (CHM). Manual measurements of crop height and biomass were taken in parallel and compared to LiDAR CHM estimates. Stepwise multiple regression was used to estimate biomass. Analysis of variance (ANOVA) tests and pairwise comparisons were used to determine di ﬀ erences between ecotypes and drought stress treatments. We found a signiﬁcant relationship between the sensor readings and manually measured crop height and biomass, with determination coe ﬃ cients of 0.73 and 0.71 for height and biomass, respectively. Di ﬀ erences in crop heights were detected more precisely from LiDAR estimates than from manual measurement. Crop biomass di ﬀ erences were also more evident in LiDAR estimates, suggesting di ﬀ erences in ecotypes’ productivity and tolerance to drought. Based on these results, application of the presented UAV-LiDAR workﬂow will provide new opportunities in assessing bioenergy crop morpho-physiological traits and in delivering improved genotypes for bioreﬁning.


Introduction
Given the rapidly growing world population and pressure placed on availability and productivity of agricultural lands by climate change, it is imperative to explore new technologies and approaches to help increase the rate of crop improvement. These improvements aim to increase biomass of food and feed crops and that of bioenergy crops which reduce the need for fossil fuels. An important step towards this goal has been the transformation in the detail and speed of genetic information brought about by next-generation sequencing. However, there remains a bottleneck in acquiring the information on the detailed morphological and physiological characteristics of large number of crop genotypes, an activity known as phenotyping. In part, this requirement has been fulfilled by the development of phenomics platforms for fast, accurate, high-resolution, and non-destructive high-throughput phenotyping [1], which have several advantages, including flexibility, convenient operation [2,3] and the ability to provide controlled experimental treatments. There have also been useful developments in the provision of high-resolution field phenotyping systems, but these are relatively large systems deployed on fixed sites [4]. However, neither of these approaches allows the flexibility to operate directly on field crops on any site required, such as a large breeding nursery, an experimental agriculture site or a typical commercial farm where detailed crop phenotyping is required. Thus, the overall aim of this study is to investigate the application of unmanned aerial vehicle (UAV)-based high-resolution remote sensing to phenotyping ecotypes of the bioenergy crop Arundo donax under contrasting water stress levels to determine the associations between genomic and phenotypic data that will allow rapid crop improvement [1,5]. Key physical parameters that are targets for A. donax high-throughput phenotypic evaluation under the real conditions experienced by the plant include crop height and above-ground biomass (AGB) production.
Formally the phenotype is an expression of the genotype and the environment in which it grows, including structural traits (e.g., plant height, leaf area index, lodging, crop canopy cover), canopy spectral texture (spectral features), physiological traits (e.g., chlorophyll, biomass, pigment content, photosynthesis), abiotic/biotic stress indicators (e.g., stomatal conductance, canopy temperature difference, leaf water potential, and senescence index), nutrients (nitrogen concentration, protein content), and yield [1]. Many of these traits can be used as indicators of plant fitness and for biomass estimation and have been used in selective plant breeding for decades. An important feature of these measurements is they often require manual measurement or sampling. This is invariably time consuming, can result in considerable crop disturbance and for some measurements requires destructive sampling. Thus, any non-invasive and non-destructive measurements which can be used as an indicator of likely biomass that can readily be applied to a biomass crop would be advantageous. Crop height is recognized as one of the most important components of AGB and yield in biomass crops [6,7]. Height can be reduced by environmental stress and can therefore be used to separate drought sensitive from drought tolerant crops [8,9]. Traditionally, crop height and diameter are measured manually with graduated poles and measuring tapes; this is difficult to do without canopy disturbance in the central region away from plot edges of tall (3.2 to 4.2 m) A. donax plants with their typically high shoot density. In terms of the choice of sensor, light detection and ranging (LiDAR) is a promising technology to measure plant height and predict biomass [10]. LiDAR equipped UAVs have been demonstrated to work well in forestry [11,12], maize [13], and shorter annual crops such as wheat [14]. There is one study using tractor mounted LiDAR to characterize a bioenergy crop, Miscanthus giganteus, with good results (reported accuracy 92-98.2%) [15]. However, the approach is limiting in terms of field size and shape as well as crop horizontal density, thus demonstrating the advantage of UAV platforms. Prior to our study, there was no survey that has applied UAV-based LiDAR to A. donax which contrasts both in terms of crop structure and leaf reflectivity with other crops studied to date.
An important aspect of LiDAR operation for agricultural and ecological studies, where timeliness of measurements is important, is the active laser illumination which means it is routinely deployed under a wide range of illumination conditions. The LiDAR data recorded is developed into a Remote Sens. 2020, 12, 3464 3 of 20 three-dimensional (3D) georeferenced point cloud further separated into a digital terrain model (DTM) and a digital surface model (DSM). These are then subtracted to create a normalized DSM, also called crop height model (CHM) when crop positions have been identified. Further analysis of the CHM allows estimation of heights and biomass [16] to answer ecologically valuable questions related to these two variables and their associated genotypes.
A. donax is a polyploid perennial rhizomatous grass grown for bioenergy production. It was found to be the most productive biomass crop when allowed to grow for more than four years [17]. Long considered an invasive plant, A. donax is unimproved and understudied as a crop, presenting a great opportunity for breeders and researchers to explore its suitability for growth on marginal lands. Given the interest in its potential for bioenergy applications, the A. donax transcriptome was sequenced and de novo assembled [18]; thus, good genetic information but much less phenotyping data is available. Not only does A. donax have high yield compared to other biomass crops, some ecotypes have high tolerance to soil contamination by heavy metals [19], tolerance to high salinity or tolerance to water stress [20,21]. Thus, it satisfies the need for an improved bioenergy crop that can yield high biomass, grow in harsh environments, and does not compete for higher quality land suitable for food production.
Here, we report the use of drone-mounted LiDAR to determine differences in responses of A. donax ecotypes to drought stress under field conditions. We set out to answer three questions: (i) can the proposed method of processing 3D point cloud data generated from LiDAR determine A. donax crop height and estimate biomass production; (ii) can it identify differences between A. donax ecotype heights and biomass; and (iii) can it identify differences between drought-stressed and well-watered ecotypes?

Plant Material and Experimental Design
The study area was located in Savigliano, northern Italy, (44 • 35 N, 07 • 37 E, 349 m above sea level) ( Figure 1a). Three clones of A. donax (named hereafter EcoA, EcoB, and EcoC) were selected from a set of 82 EuroMediterranean ecotypes [22] and vegetatively propagated by rhizomes. The selected ecotypes originated from three contrasting environments: EcoA was collected from the coastal habitat of Attica region (southern Greece), characterized by typical warm Mediterranean climate; EcoB originated from the coastal habitat of southern Dalmatia region (southern Croatia), characterized by temperate oceanic climate; and EcoC was collected from the hilly area of Viseu district (northern Portugal), characterized by temperate Mediterranean climate. The field experiment began in March 2014. There were three replicate plots of each ecotype grown in each of two watering regimes: the well-watered (WW, kept at 40% volumetric water content) and natural moderate drought stress (mDr, allowed to fall to 15% volumetric water content during June). Each ecotypic replicate consisted of 30 homogeneous cuttings which were planted at distances of 0.5 m × 1.0 m. The experimental design is shown in Figure 1.

UAV Phenotyping Platform
In our study, the UAV LiDAR data were acquired by an unmanned G4 SkycraneV2 ® octocopter (Multirotor Service-drone, Germany) equipped with micro LiDAR Velodyne HDL-32E microsensor (Velodyne, San Jose, CA, USA) ( Figure 2). The hyperspectral sensor was not used in this study. The micro LiDAR sensor is small and lightweight, featuring up to ±2 cm accuracy, 40 • vertical field of view and 360 • horizontal field of view with 32 laser beams/scan and dual returns up to 70 m (Velodyne Acoustics, Inc., Morgan Hill, CA, USA). An inertial measurement unit is combined with a real time kinematic global positioning system (GPS) to help georeferencing the captured images, with a precision of 1.0 and 2.0 m in the horizontal and vertical directions, respectively.
Remote Sens. 2020, xx, x 5 of 20 In our study, the UAV LiDAR data were acquired by an unmanned G4 SkycraneV2 ® octocopter (Multirotor Service-drone, Germany) equipped with micro LiDAR Velodyne HDL-32E microsensor (Velodyne, San Jose, CA, USA) ( Figure 2). The hyperspectral sensor was not used in this study. The micro LiDAR sensor is small and lightweight, featuring up to ± 2 cm accuracy, 40° vertical field of view and 360° horizontal field of view with 32 laser beams/scan and dual returns up to 70 m (Velodyne Acoustics, Inc., Morgan Hill, CA). An inertial measurement unit is combined with a real time kinematic global positioning system (GPS) to help georeferencing the captured images, with a precision of 1.0 and 2.0 m in the horizontal and vertical directions, respectively.

UAV-LiDAR Campaign
The airborne campaigns were conducted on 29 June 2016, with UAV LiDAR scanning at an altitude of 25 m above ground with separate flight plans for each treatment (Figure 1a). Four ground control points were installed for each treatment (Figure 1a), and their positions established with an onboard GPS and an on-ground global navigation satellite system with an accuracy of 3 mm (Leica GS08plus, Leica Geosystems, Switzerland) operating in real-time kinematics.
The flight mission was planned using the autopilot Ground Station v.3.2.304 (Multirotor Service-Drone, Germany). The drone was flown autonomously (GPS-waypoint navigation) at a nominal speed of 3 m/s. Experimental plots were scanned for two 9-min flights under stable cloudless and low-wind conditions. Flights were performed at 14:22 local time above the WW and at 15:29 local time above the mDr field. In order to capture each experimental plot in a single flight pass, the flight plan consisted of transects parallel to the plant rows ( Figure 1a). A ground station processed the UAV safety manual control and sent telemetry data (position, attitude, and status data) through a radio link at 2.4 GHz to a field laptop. Figure 3 shows a UAV RGB image of A. donax ecotype in the field and 3D LiDAR point clouds.

UAV-LiDAR Campaign
The airborne campaigns were conducted on 29 June 2016, with UAV LiDAR scanning at an altitude of 25 m above ground with separate flight plans for each treatment (Figure 1a). Four ground control points were installed for each treatment (Figure 1a), and their positions established with an onboard GPS and an on-ground global navigation satellite system with an accuracy of 3 mm (Leica GS08plus, Leica Geosystems, Switzerland) operating in real-time kinematics.
The flight mission was planned using the autopilot Ground Station v.3.2.304 (Multirotor Service-Drone, Germany). The drone was flown autonomously (GPS-waypoint navigation) at a nominal speed of 3 m/s. Experimental plots were scanned for two 9-min flights under stable cloudless and low-wind conditions. Flights were performed at 14:22 local time above the WW and at 15:29 local time above the mDr field. In order to capture each experimental plot in a single flight pass, the flight plan consisted of transects parallel to the plant rows ( Figure 1a). A ground station processed the UAV safety manual control and sent telemetry data (position, attitude, and status data) through a radio link at 2.4 GHz to a field laptop. Figure 3 shows a UAV RGB image of A. donax ecotype in the field and 3D LiDAR point clouds.
low-wind conditions. Flights were performed at 14:22 local time above the WW and at 15:29 local time above the mDr field. In order to capture each experimental plot in a single flight pass, the flight plan consisted of transects parallel to the plant rows ( Figure 1a). A ground station processed the UAV safety manual control and sent telemetry data (position, attitude, and status data) through a radio link at 2.4 GHz to a field laptop. Figure 3 shows a UAV RGB image of A. donax ecotype in the field and 3D LiDAR point clouds.

Plant Height and Biomass Ground Truthing
The heights of four random stems (StH [cm]), in the central 1 m 2 of each plot, were manually measured to assess the significance of their relationship with the UAV-LiDAR data at the same time the UAV was flying, forming the ground-truth dataset. Heights were measured from the basal node to the top node with a graduated telescopic pole. Nodes are distinct structural components of the stem where the leaves grow out. They are easy to be recognized at a distance by the naked eye as there will be a tiny knob-like on the stem. The lowest node, which is usually just above the ground in A. donax is used as the lower reference point for our stem height measurements. For the most part of the A. donax growth, the position of the uppermost node is hidden because of overlapping leaf bases and used for the stem height measurement to avoid human bias in the selection of the tallest leaf at the top which is much harder to determine.
Finally, at the end of the growing season (November 2016), the time when A. donax is usually harvested as a bioenergy crop, above-ground biomass (total dry mass, AGB, g m −2 ) was estimated for all replicates. Stems were cut at 5 cm above ground level and weighed to determine fresh weight. Then, sub-samples of mixed tissues (stems, leaves, and flowers) were oven-dried at 75 • C for 72 h to determine the conversion to dry weight.

Data Processing and CHM
The CHM from the LiDAR data was created according to the workflow in Figure 4. The first step was the conversion of each UAV LiDAR scan from .pcap to .las format files using LidarTools (Hyperspec III, version 3.1, Headwall Photonics Inc., Boston, MA, USA). The second step was to clip the LiDAR point clouds to the study areas using LASclip algorithm in LAStools and then eliminate noise (first and last pulse) and outliers inside these point clouds using the LASnoise algorithm in LAStools [23]. LASground, from LAStools, was then used in the third step, to separate ground returns from non-ground returns. DTMs and DSMs with grid cell size of 20 cm for each field treatment were generated using the "lidR" R package [24]. In the fourth step, DTMs and DSMs were processed in the geographic information system application QGIS (QGIS Development Team, 2018) by normalizing the DSM with the DTM to obtain the CHM (that is the residual distance between the ground and the top of the crop above the ground or the DSM minus the DTM). The CHM was georeferenced by way of ground control points. Finally, an area-based approach [16] was used to generate height and biomass prediction in an area of 1 m 2 around the center of each ecotypic replicate as described in detail below.

CHM Validation and AGB Estimation
LiDAR heights were extracted from a 1  1 m region at the center of each replicate ( Figure 3a) and averaged for validation against the four randomly measured stems within that central square (Figure 1c). The relationship between LiDAR-estimated and manually measured crop heights was evaluated using a linear regression analysis. The coefficient of determination (R 2 ) was used to define the model's predictive accuracy (considering an F-test p-value < 0.05 as significant and p-value < 0.001 as highly significant), and the root mean square error (RMSE), as a metric for accuracy. Statistical analyses were conducted using R packages such as the "leaps" and the "metrics". While the leaps performs an exhaustive search for the best subsets of the variables in x for predicting y in linear

CHM Validation and AGB Estimation
LiDAR heights were extracted from a 1 × 1 m region at the center of each replicate ( Figure 3a) and averaged for validation against the four randomly measured stems within that central square (Figure 1c). The relationship between LiDAR-estimated and manually measured crop heights was evaluated using a linear regression analysis. The coefficient of determination (R 2 ) was used to define the model's predictive accuracy (considering an F-test p-value < 0.05 as significant and p-value < 0.001 as highly significant), and the root mean square error (RMSE), as a metric for accuracy. Statistical analyses were conducted using R packages such as the "leaps" and the "metrics". While the leaps performs an exhaustive search for the best subsets of the variables in x for predicting y in linear regression [25], the metrics implements metrics for regression, time series, binary classification, multiclass classification, and information retrieval problems [26].
For the AGB estimation, eleven LiDAR CHM metrics (defined in Table 1), calculated for 1 m 2 area of each plot, were assessed as potential variables for inclusion in the model. Following the area-based approach [16], these metrics were computed for the 1 m 2 area of each ecotypic replicate. We used several linear regressions to study the relationship between each CHM metric variable, extracted from LiDAR data, and the manually measured AGB. Initially, LiDAR CHM metrics with Pearson correlation coefficient (R) values larger than 0.4 were selected to estimate AGB, as is usually done to reduce the number of possible AGB estimators and speed up computations (e.g., [27][28][29]). The correlation coefficients of all the LiDAR metrics and the field measured AGB are plotted in Figure S1. Two way-ANOVA tests and pairwise comparison with Bonferroni correction (to adjust p-values for multiple comparisons) [30] were performed on the ground data and the LiDAR estimates to determine whether differences between ecotypes and water treatments exist. Linear regressions between AGB and heights were also evaluated for manual measurements and LiDAR estimates.

LiDAR CHM
The LiDAR system generated high-density 3D data in discrete point clouds ( In this study, there was a highly significant positive linear relationship (R 2 = 0.73, p-value < 0.001) between the estimated A. donax crop heights (LiDAR Hmean) derived from LiDAR CHM ( Figure 6) and the manually measured crop heights (Figure 7a). No significant difference was found between the relationships of the two treatments ( Figure 7b) indicating that the relationship between LiDAR estimated crop heights and ground measurement is stable. Although the relationships between Hmean and manual height measurements were significant, the Hmean values were consistently greater. This could be related to the greater number of points identified as stems (12-24; Figure S1) compared with the four stems measured manually. However, difference in heights is more likely to result from different detection of the top of the stem between manual and LiDAR measurements. It can be seen in Figures 3 and 6 that the LiDAR returns are apparently from upper leaves in the canopy, which were not included in manual stem height measurements as mentioned in the methods Section 2.3.2. In this study, there was a highly significant positive linear relationship (R 2 = 0.73, p-value < 0.001) between the estimated A. donax crop heights (LiDAR Hmean) derived from LiDAR CHM ( Figure 6) and the manually measured crop heights (Figure 7a). No significant difference was found between the relationships of the two treatments ( Figure 7b) indicating that the relationship between LiDAR estimated crop heights and ground measurement is stable. Although the relationships between Hmean and manual height measurements were significant, the Hmean values were consistently greater. This could be related to the greater number of points identified as stems (12-24; Figure S1) compared with the four stems measured manually. However, difference in heights is more likely to result from different detection of the top of the stem between manual and LiDAR measurements. It can be seen in Figures 3 and 6 that the LiDAR returns are apparently from upper leaves in the canopy, which were not included in manual stem height measurements as mentioned in the methods Section 2.3.2. Remote Sens. 2020, xx, x 10 of 20

AGB Estimates Derived from 3D LiDAR Point Clouds
None of the individual LiDAR-derived metrics had a highly significant relationship with AGB (Table 2) with the highest R 2 being 0.34. However, six of the eleven metrics showed a small but significant (R 2  0.22, p-value < 0.05) relationship with manually measured AGB (Table 2). Hsum, Hmin, Hrange, Hmajority, and Hvariety did not have a significant relationship (p-value > 0.05) with AGB.
A stepwise multiple linear regression was then used to eliminate multicollinearity issues and to select the optimal regression equation. The LiDAR metrics, Hmax, Hmajority, and Hvariety, gave the best model to predict total biomass with R 2 of 0.71 and RMSE of 908.5 g m −2 (Figure 8). The optimal regression equation is defined as: where AGB is above-ground biomass (g m −2 ), Hmax is the maximum of stem height values in cm, Hmajority is the stem height with most occurrences in cm, and Hvariety is the count of unique stem height values.

AGB Estimates Derived from 3D LiDAR Point Clouds
None of the individual LiDAR-derived metrics had a highly significant relationship with AGB (Table 2) with the highest R 2 being 0.34. However, six of the eleven metrics showed a small but significant (R 2 ≥ 0.22, p-value < 0.05) relationship with manually measured AGB (Table 2). Hsum, Hmin, Hrange, Hmajority, and Hvariety did not have a significant relationship (p-value > 0.05) with AGB.
A stepwise multiple linear regression was then used to eliminate multicollinearity issues and to select the optimal regression equation. The LiDAR metrics, Hmax, Hmajority, and Hvariety, gave the best model to predict total biomass with R 2 of 0.71 and RMSE of 908.5 g m −2 (Figure 8). The optimal regression equation is defined as: where AGB is above-ground biomass (g m −2 ), Hmax is the maximum of stem height values in cm, Hmajority is the stem height with most occurrences in cm, and Hvariety is the count of unique stem height values.

Heights and Biomass of Three A. donax Crop Ecotypes Under Natural Moderate Drought
The crop height manual measurement in WW treatment ranged from 226.5 cm for EcoB to 369.3 cm in EcoA with a mean of 307.5 cm, while the mDr treatment varied from 237.8 cm in EcoB to 350 cm in EcoA with a mean of 298.2 cm (Figure 9a,b). ANOVA test showed that the differences between the ecotypes were significant, F(2,12) = 24.47, p-value < 0.001. No differences were observed across treatments, with F(1,12) = 0.38, p-value = 0.54; the interaction of ecotype with water treatment was also not significant, F(2,12) = 0.48, p-value = 0.63. Doing a pairwise comparison we found that EcoB is significantly shorter than the other two ecotypes, p-value < 0.001 (Figure 9a). For the LiDAR estimated

Heights and Biomass of Three A. donax Crop Ecotypes under Natural Moderate Drought
The crop height manual measurement in WW treatment ranged from 226.5 cm for EcoB to 369.3 cm in EcoA with a mean of 307.5 cm, while the mDr treatment varied from 237.8 cm in EcoB to 350 cm in EcoA with a mean of 298.2 cm (Figure 9a,b). ANOVA test showed that the differences between the ecotypes were significant, F(2,12) = 24.47, p-value < 0.001. No differences were observed across treatments, with F(1,12) = 0.38, p-value = 0.54; the interaction of ecotype with water treatment was also not significant, F(2,12) = 0.48, p-value = 0.63. Doing a pairwise comparison we found that EcoB is significantly shorter than the other two ecotypes, p-value < 0.001 (Figure 9a). For the LiDAR estimated crop height, we also found similar results, but in addition, doing the pairwise comparisons we found a significant effect of drought treatment on the height of EcoA (p-value < 0.05).
Remote Sens. 2020, xx, x 13 of 20 treatments p-value < 0.05 (Figure 9c). Finally, with regards to LiDAR estimated biomass, we found significant effect of ecotype, drought treatment, and their interaction; F(2,12) = 4.43, p-value<0.05, F(1,12) = 17.28, p-value = 0.001, F(2,12) = 5.61, p-value < 0.05, respectively (Figure 9d). EcoA was again much more productive in the WW plot. Ecotypes situated in WW were on average significantly more productive than ecotypes in mDr. EcoB was found to be less productive than the other two ecotypes; however, because of the interaction between ecotype and treatment, EcoB was not showing any sign of stress in mDr, while EcoC did show some insignificant stress effect (Figure 9d). Looking at crop biomass and crop height, there was a positive but non-significant relationship in the manually measured dataset (Figure 9e). This relationship was found to be significant when looking at LiDAR data (Figure 9f).  The AGB, obtained from the harvesting of stems in each treatment, ranged from a minimum of 2096.3 g m −2 (EcoB in WW) to a maximum of 7845.9 g m −2 (EcoC in WW) with an average of 4837.6 g m −2 . ANOVA tests showed mainly no significance between ecotypes, water availability treatments and their interaction. We did find, however, a significant difference in EcoA AGB between the treatments p-value < 0.05 (Figure 9c). Finally, with regards to LiDAR estimated biomass, we found significant effect of ecotype, drought treatment, and their interaction; F(2,12) = 4.43, p-value<0.05, F(1,12) = 17.28, p-value = 0.001, F(2,12) = 5.61, p-value < 0.05, respectively (Figure 9d). EcoA was again much more productive in the WW plot. Ecotypes situated in WW were on average significantly more productive than ecotypes in mDr. EcoB was found to be less productive than the other two ecotypes; however, because of the interaction between ecotype and treatment, EcoB was not showing any sign of stress in mDr, while EcoC did show some insignificant stress effect (Figure 9d).
Looking at crop biomass and crop height, there was a positive but non-significant relationship in the manually measured dataset (Figure 9e). This relationship was found to be significant when looking at LiDAR data ( Figure 9f).

Discussion
To date, LiDAR-derived metrics have been used to predict attributes such as AGB in forests and crops [7,14,16,[31][32][33][34][35]. Given the lack of fast and non-destructive alternatives, the capacity to estimate AGB using UAV-derived LiDAR metrics is a critical outcome for dense bioenergy crops such as A. donax. Whilst physical sampling is the most accurate method for estimation of AGB, it is invasive, destructive, costly, and can only be used to study the time course of AGB development in large plots with extensive sampling. The UAV-mounted LiDAR sensor was developed for large field experiments and breeding trials and its deployment has clear advantages over current practice [2]. It is easily transported to any required location overcoming a major limitation of fixed phenotyping platforms. It is significantly more flexible than using airplanes where experiments are constrained by location, cloudiness, and expense [4,36]. This study has demonstrated for the first time the practicality of this approach for the bioenergy crop, A. donax.

Crop Height Estimation
Crop heights derived from the LiDAR measurements in A. donax had an R 2 of 0.73 with manual height measurements (Figure 7). This result is similar to R 2 for other crop height measurements estimated from LiDAR, of sugar beet (R 2 = 0.70), wheat (R 2 = 0.78), potato (R 2 = 0.50) and sorghum (R 2 = 0.63) [14]. Although these results demonstrate that LiDAR-derived height can estimate the differences in crop heights across ecotypes, a systematic difference between the manual stem height measurement and the UAV-measured height was observed (Figure 7). One possible explanation for the differing results is that the complex canopy surface makes it difficult to compute a high-quality point cloud [37]. Moreover, in practice, most LiDAR returns are likely to come from leaves rather than the denser structures containing more of the biomass such as the stems. As previously mentioned, measuring A. donax plant height on the ground by eye is challenging because of the density and height of the stems, and hence to avoid bias in measurement, we standardized the process by measuring each stem from the basal node to the top node. Another explanation is that the LiDAR ground penetration is hindered by the thickness of the canopy [35,38] or by the footprint of the LiDAR itself.

AGB Estimation
Linear relations between individual LiDAR-derived-CHM metrics and biomass had low determination coefficients and high mean square errors ( Table 2). As previously mentioned, these might be related to the poor penetration of the laser within the canopy [35,38]. However, an individual metric is unlikely to adequately represent the complex canopy structure. There are also other reports of low correlations between LiDAR and manual measurements [39]. Such limitation could be reduced by using a LiDAR with a smaller footprint to better penetrate the leafy biomass. Another way to improve biomass prediction is to take into consideration several other crop cover metrics, which may represent other features of canopy structure. Here, from the eleven LiDAR-derived metrics, Hmax, Hmajority, and Hvariety were together the best variables to include in a multiple regression for biomass estimation (Equation (1)) with good coefficient of determination (R 2 = 0.71) and low bias (bias = 3.2) (Figure 8). In the 1 m 2 evaluated these variables represent respectively the maximum height captured, the mode of the points and the count of distinct values. As these variables greatly improved the estimate, they may be recovering other features that contribute to canopy volume. While Hmax contributed the most to prediction of biomass, Hmajority and Hvariety individually correlated negatively with the dry biomass estimation, meaning that tall stems contribute more to the yield than the height with the most occurrences in the crop, and that crops which distribute resources more uniformly across stems (smaller Hvariety) yield more biomass. The positive linear relationship obtained is higher than that found for sugar beet (R 2 = 0.68), lower than that for the much denser crop of winter wheat (R 2 = 0.82), better than that for potato (R 2 = 0.24) [14], and similar to that obtained for different nitrogen treatments in barley (R 2 = 0.71) [40].
In our study, A. donax biomass was measured in November, compared with the earlier acquisition of LiDAR data in late June. Since stem height mostly contributes to AGB, several studies have shown that A. donax reaches the height plateau around July-August [41][42][43]. In addition, to maximize biomass yield, A. donax harvesting in North Italy is usually done in November just before the crop starts to senesce [22,44]. Waiting this long to collect LiDAR data is however not an option, as the inflorescences produced in August can significantly affect crop height and make it difficult to acquire accurate data manually. Despite the difference in timing between collection of height data and harvesting, multivariate LiDAR-derived biomass estimates related well to the dry weight ( Figure 8). Biologically this does not come as a surprise, as stems produced in June are the most productive in terms of yield and later stems produced are smaller. If this relationship proves robust and can be established before flowering, it may have useful predictive ability for selecting suitable material for crossing in breeding programs.
Finally, the relationship between crop biomass and height was significantly better when looking at LiDAR estimates than at manual measurements (Figure 9e,f). The difference between the measured biomass and the estimated biomass can be explained by number of stems weighed (four stems) as opposed to the number of stems present in the 1 m 2 detected by the LiDAR. The difference between biomass-height relationship, on the other hand, is probably due to the inaccuracy of human judgment in measuring crop aspects as opposed to LiDAR biomass being derived from the LiDAR crop heights. This is a specific problem in this very tall (3.2 to 4.2 m) dense crop where it is difficult to see the top of the stems, making it really hard to measure their heights without introducing errors. A greater number of plants should be surveyed by LiDAR with accompanying ground truth data in future studies to improve the crop biomass-height relationship and determine whether it is linear or exponential for A. donax.

Comparison of Ecotypes
These data confirm that EcoB, originally from a temperate oceanic climate, was significantly shorter than the other two Mediterranean ecotypes of A. donax (Figure 7a). Significant differences in height and biomass were detected in EcoA between WW and mDr from LiDAR estimates (Figure 9b). While ground height measurement showed no effect of drought stress on EcoA (Figure 9a), manually measured biomass did (Figure 9c). This difference in EcoA response to drought stress was also much more significant when looking at LiDAR estimated biomass (Figure 9b,d). This suggests that LiDAR phenotyping was better at determining differences in biomass and height than manual measurements, and LiDAR is currently the only non-destructive alternative. The much higher spatial sampling possible in the LiDAR derived CHM compared with the manual sampling makes LiDAR the preferred approach. Similar height measurement inaccuracy was reported in sorghum [45], and it suggests that it is difficult without a very large investment in time and care to record accurately the distribution of height in these tall field crops.
The lack of differences in EcoB and EcoC between treatments suggests that these A. donax ecotypes are more tolerant than EcoA to moderate drought stress applied during the experiment. There are other reports of stress tolerance in A. donax [19][20][21]46], hence its suitability for growth on marginal lands. Future work will aim to increase drought intensity and duration to determine how much is A. donax growth reduced at different drought stress levels and which ecotypes are more tolerant to these stress levels.
Our findings demonstrate that UAV-LiDAR is able to capture phenotypic differences across the three A. donax ecotypes, and therefore can be applied to larger scale than currently performed experimental trials monitoring the productivity of this bioenergy crop under different environmental conditions (e.g., [47]). The differences in productivity suggest that EcoA is much more productive than EcoB and EcoC. However, given EcoA's low tolerance to drought, it is clear that it might not be the most suitable for bioenergy crop production on marginal lands. EcoB, on the other hand, had the lowest biomass but showed no difference in performance between the two treatments, which suggest that EcoB is more adapted to water scarcity in drought-prone environments, can grow on marginal lands and might yield more biomass than EcoA under similar soil conditions. EcoC showed intermediate performance for biomass production and drought tolerance, indicating a potential for bioenergy crop production on marginal lands. More research on these ecotypes is needed to gain a deeper understanding of the contributions of phenology, physiology, morphology to biomass production in this important bioenergy crop. However, depending on whether we want to maximize biomass yield or drought tolerance, the ecotypes studied demonstrate a huge potential for each scenario.

Implications
Future work will use more ecotypes at several locations and time points to further verify the robustness of the positive linear relationship between LiDAR estimates and biomass or other biophysical traits at various growth stages in response to environmental stresses. Sensors providing hyperspectral and thermal data, for example, will also be tested for their ability to contribute to novel associations between the phenotype and genotype of A. donax that are of interest to breeders and growers. Monitoring the productivity during the growing season can also be very important to optimize the biomass yield. Although the results in crop height and biomass estimation have demonstrated the potential of the UAV-LiDAR platform, some limitations are noted and should be worked on to increase utility of the phenotyping platform. For instance, UAV-LiDAR derived height measurements systematically overestimated the height of the A. donax. This is acceptable for comparison of crop heights as used in many breeding experiments; however, for absolute height measurements and more accurate biomass estimation, developing robust calibration with ground data should be the eventual aim.

Conclusions
Recent efforts in A. donax crop improvement are focused on the measurement of AGB and crop height as determinants of greater biomass yield. We demonstrated that the present methodology with the UAV-mounted LiDAR sensor is particularly useful for high-throughput phenotyping of A. donax in the field; it has a high potential for advancing A. donax crop phenomics and genomics-assisted breeding.
Our findings could help biomass crop breeders and researchers to reduce the gap between phenomics and genomics by generating accurate A. donax information in an efficient, non-destructive, and inexpensive way when compared to traditional laborious approaches. The methodology described also has potential for data mining to select the most promising ecotypes, making it possible to assess specific multiple criteria for ranking ecotypes, one of the main techniques used in crop breeding. This would be a valuable extension of the crop height and AGB successfully estimated from LiDAR in A. donax. We were able to capture differences in crop heights between treatments and between ecotypes that manual measurements were unable to detect. This improvement was probably related to the unbiased sampling at 30 or more points per m 2 of the crop compared with the four stems selected and measured manually. These results are the first demonstration of the suitability of this UAV-based workflow for A. donax field phenotyping. Our findings show that UAV-based LiDAR scanning can replace traditional field-based crop surveys of mean height and AGB of different ecotypes and under different environmental conditions. This methodological and technological tool could be adapted to provide technical support to promote commercially viable applications of UAVs in crop phenotyping of other bioenergy crops (e.g., miscanthus and switchgrass) for other research programs and breeding companies.
Future investigations could be focused on the development of comparative analyses of A. donax ecotypes to test their growth habits and AGB yields over multiple environments, locations, and years to select the best performing and stable ecotypes. Considering that our research study was on three A. donax ecotypes stressed at one time point, next work could focus on time-series measurements of a large ecotypic panel, since breeders may need to study a large pool of plant materials and produce ecotypes that are also suitable for use in multiple cropping systems.
Supplementary Materials: The following will be available online at http://www.mdpi.com/2072-4292/12/20/3464/ s1, Figure S1. Correlation matrix plot with significance levels between the different light detection and ranging (LiDAR) metrics and the field measured above-ground biomass (AGB) production in Arundo donax. The lower triangular matrix is composed of the bivariate scatter plots with fitted smooth red lines. The upper triangular matrix shows the Pearson correlation significance level (as red asterisks or red squares). Each significance level is associated to a symbol: p-values 0.001 (***), 0.05 (*), 0.1 ( ).